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ABSTRACT 


Geoid undulation computations were carried out at 39 laser stations 
distributed around the world using a combination of terrestrial gravity 
data within a cap of radius 2 ' and a potential coefficient set (Rapp and 
Cruz, 1986b) up to degree 180. The traditional methods of Stokes’ and 
Meissl’s modification together with the new Molodenskii’s method and the 
modified Sjoberg's method were applied. Performing numerical tests 
based on global error assumptions regarding the terrestrial data and 
the geopotential set we concluded that the modified Sjoberg’s method is 
the most accurate and promising technique for geoid undulation 
computations. The numerical computations of the geoid undulations 
using all the four methods resulted in agreement with the "ellipsoidal 
minus orthometric" value of the undulations on the order of 60 cm or 
better for most of the laser stations in the eastern United States, 
Australia, Japan, Bermuda and Europe. A systematic discrepancy of 
about 2 meters for most of the western United States stations was 
detected and verified by using two relatively independent data sets. 
The cause of this discrepancy was not found. A correction due to the 
inconsistencies of the terrestrial data and the potential coefficients 
within the cap surrounding the laser station, called the "local average" 
correction improved the results by 30 cm and it seems necessary to 
apply this correction. For oceanic laser stations in the western Atlantic 
and Pacific oceans that no terrestrial data available, the adjusted 
GEOS-3 and SEASAT altimeter data were used for the computation of the 
geoid undulation in a collocation method. 
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CHAPTER I 


INTRODUCTION 

The purpose of this investigation is the precise absolute geoid 
undulation computation for certain laser stations distributed around the 
world. The classical methods that combine terrestrial gravity data 
within a cap surrounding the computation point together with a high 
degree potential coefficient set through the use of integral formulas 
have been implemented to attack the problem. The only exception is the 
computation of the geoid undulation at four laser stations in the Pacific 
and two stations in West Atlantic that have been carried out using 
adjusted Geos-3 and SEASAT altimeter data in a collocation method. 
Chapter 2 will give the necessary truncation theory. The general 
formulation is given in section 2.1. The conventional method of Stokes’ 
and Meissl’s modification are presented in Section 2.1 as a review. The 
modification of Stokes' function with the ultimate goal to minimize certain 
errors in geoid undulation computations is a common concept for both 
Molodenskii’s and modified Sjoberg's methods. This is why Molodenskii’s 
and modified Sjoberg’s methods are discussed together as "new" 
methods in Section 2.3. 

After the theoretical background for the above four methods is 
established, numerical computations based on certain error models will 
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show us the behavior of the global root mean square (RMS) errors of 
the four methods as a function of the capsize. This is discussed in 
Chapter 3. We then will present information in Chapter 4 concerning 
the selected laser stations around the world. The available gravity data 
surrounding the laser stations were processed to be used in the 
integration formulas. This process is given in Chapter 5. The 
pre-processing of the data is first discussed in Section 5.1. Some 
prediction techniques which are necessary to transform the point data 
to mean values will be given in Section 5.2. The prediction procedure 
we used for the Bermuda area will be separately described in Section 
5.3. The terrain corrections and the indirect effect for stations on high 
mountains is discussed in Section 5.4. To avoid the spherical 
approximation, ellipsoidal correction formulas will be given for the 
conventional methods, and will be derived for the new methods, in 
Chapter 6. A local average correction which is new in concept will be 
given in Chapter 7. 

The numerical results for geoid unduation computations using all the 
four methods (Stokes’, Meissl’s, Molodenskii’s and modified Sjoberg’s 
methods) will be presented in Chapter 8: The details for the numerical 

application of the four methods will be given in Section 8.1; the results 
will be given in Section 8.2; a discussion of the results and comparisons 
will be given in Section 8.3. Emphasis will be given on the following 
comparisons: 

a) Terrain corrected vs. uncorrected undulations. 

b) Computed undulations using all the four method vs. the value of 
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the undulations computed as the difference ellipsoidal minus 
orthometric height. 

c) Undulations computed using two different methods. 

d) Undulations computed using the four methods and the values of 
the undulations taken from a high degree geopotential model. 

e) Undulations computed excluding the local average correction vs. 
undulations computed including the local average correction. 

The results will be presented by regions, for detection of any 
systematic differences that are correlated with a specific region. 
The collocation method was also used to compute undulations from 
altimeter data in the vicinity of Bermuda. These results together with 
the undulations obtained from the calibration of SEASAT will be 
compared to the gravimetric undulations and will be presented in 
Section 8.4. The undulations for five oceanic laser stations computed 
using altimeter data in a collocation method will be given in Section 8.5. 
Finally, summary and conslusions will be given in Chapter 9, and 
recommendations for future studies in Chapter 10. 

An attempt has been made to maintain a uniform simple notation. 
The following will be adhered to: 

N geoidal undulation 
R mean earth radius 

7 mean value of the normal gravity 

Y>o spherical cap of radius Vo surrounding the computation 


4 


S(cos^) 

S 0 

W, (cosi/>) 
i=l 
i=2 
i=3 
i=4 
S i (cos^) 
M 
n 

n T 


hi 


4gg 


M(-) 

C n =M[ ( Ag n ) 2 ] 

<5C n =M[ (* = )*] 
*n=M[U!) 2 ] 

Qi n 


point 

Stokes’ function 
S ( cosV'o ) 

proper function subtracted from S(cos V) 

Stokes ’ method 
Meissl’s method 
Molodenskii’s method 
Modified Sjoberg method 
S(cosV) - Wj (cos^O 

maximum degree of potential coefficients used 
maximum number of Fourier coefficinets of W f (cosY>) 
maximum degree of the terrestrial error degree 
variances 

n th surface spherical harmonic of terrestrial Aj given 
on a sphere of radius R 

n th surface spherical harmonic of A|, based on a set of 
known potential coefficients given on a sphere of 
radius R 

degree error of Aj n (A^n = 4^ + ej) 
degree error of A® (Ag® = A| n + e®) 
global average operator 

degree variance of Ag n for high degrees given on a 

sphere of radius R 

degree variance of 

degree variance of eJ 

truncation coefficients 



X in Fourier coefficients of S ( (cosy) 

W, n Fourier coefficients of W^cosf) 

Legendre polynomial of degree n 


P„(cosV) 



CHAPTER II 


TRUNCATION THEORY 


2.1 General Formulation 

The geoid undulation N with respect to a specified mean earth 
ellipsoid can be computed using gravity anomalies Ag on the geoid 
through the use of the Stokes* integral: 

N = ~ JJ S(cos^)Agd<r (2.1) 

a 

A number of assumptions are associated with formula (2.1): 

a) The boundary surface of the geoid is a sphere <r 

b) Integration is done using global gravity data (on the sphere a) 

c) No masses external to the geoid exist 

d) The mass of the reference ellipsoid equals the earth's mass 

e) The normal potential on the ellipsoid equals the gravity potential 
on the geoid 

f) The ellipsoid's center coincides wtih the earth's center of mass. 

Assumption a) can be avoided by computing the ellipsoidal 
corrections (see Chapter 6); the error when doing assumption c) can be 
diminished by the use of the terrain-corrected gravity anomalies (see 
Section 5.4); assumptions d), e) and f) cannot be avoided and if they do 
not hold, a zero and first order term will have to be added to the 
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computed undulation. Finally, assumption b) can be more or less 
avoided if gravity anomalies within a cap <r c surrounding the 
computation point and a given geopotential model are combined. All four 
methods that are described below will split the integral ( 2 . 1 ) into two 
parts (gravity anomalies within a cap and potential coefficient 
information) to avoid assumption b) which cannot be rigorously avoided 
due to the lack of global coverage of the existing gravity data. 

The next Btep is to consider the basic modification of (2.1), so that 
certain errors are minimized. For all the four methods, denoting with 
S | (cos^) the modified Stokes’ function, the modification is done so that: 
(i = 1, 2, 3, 4 corresponds to the four methods: Stokes’, Meissl’s, 

Molodenskii’s, modified Sjoberg’s). 

St(cosY0 = S(cosV') - Wj (cos Y>) (2.2) 

where W^cos^ 1 ) is a properly defined function of the spherical distance 
i>. [For the original Stokes’ method it is clear that no modification is 
done, but the notation in ( 2 . 2 ) will be kept even for i = 1 , for 
convenience. It simply holds then W x (cos i>) - 0; Si(cosV) = S(cosY')]. 

Any modification of the Stokes’ function has to maintain the 
integrity of the basic Stokes’ equation in (2.1). Substituting (2.2) into 
(2.1), the modified Stokes’ equation will be: 

N< = 4^ JJ S, (cosf)Agdff + 4 ^; JJ W, (cosY-)Agd<r (2.3) 

<r (r 

Equation (2.3) is rigorously equivalent to (2.1) 
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A powerful relationship that provides an immediate conversion of an 
integral over the sphere a to a series of harmonics is proven in (Jekeli, 
1980) and will be used extensively within the text. The relationship 
reads: 

JJ K(cosi»Ag nm d a = 2wk n Ag nm (2.4) 

a 

where K(cosY') an arbitrary kernel function 

27rk n the eigenvalues of the integral operator in (2.4) 

Ag nffl the n th degree and m th order spherical harmonic of the 
surface function Ag. 

From (2.4) we obtain 

X I IT K(cosV)Ag nin d<r = X I 27rk n Ag nm , or 

n-o m=o n-o m-o 

K(cos^)Agdff = 2n Z k n Ag n (2.5) 

n=o 

Applying (2.5) to (2.3) we easily obtain: 

N, = 4^ JJ s l(cosV)Agd<r + ^ WlnAgn ( 2.6) 

<f 

where 2nVf t n are the eigenvalues of the second integral operator in 
(2.3). 



Now we are ready to split (2.6) (or equivalently (2.1)) into two 
components. From (2.6) we take 


N = 4^ JJ S,(cosf)Agd<r + jj S< (cosVOAgdff + ^ J 0 W in*gn 

c G'-G c 


(2.7) 
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Note that the second and the third term in (2.7) are to be computed 
from a given set of potential coefficients and the first term is the cap 
contribution. 


To bring the second integral to the form of a summation, we write 


4 ~ II s i (cos^O Agd<r = JJ S < (cos?) Agd<r 

l CT-<T i ' <r 


( 2 . 8 ) 


where McosY) = { 0 if 0 ‘ * ‘ ♦« 


S f (cos?) if ? c < V * w 


(2.9) 


Expanding S^cost/i) in spherical harmonics we have 


2n+l 


Sj(cOS?) = l o - Qi nP n (COS^) 

n = o ^ 


( 2 . 10 ) 


where Q 1n are the Fourier coefficients of S j (cosV'): 


Qln = 2^ II s i (cos?)P n (cos?)da 


( 2 . 11 ) 


Since 2wQ ( n are the eigenvalues of the integral operator of the 
right-hand side of equation (2.8), equation (2.5) together with (2.8) 
gives 


4 “ II Sf (cos?)Agd<r = ^ J o Q jn Ag n (2.12) 

Substituting (2.12) into (2.7) the final equation for the modified Stokes’ 
function is 

N| = if; II S,(cosY>A*dw + OlA*. + jS £ W,„Ag„ 

G r 
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= 4 ^ If S,(cos^)Agd ff + ^ i(Q 1n + W in )Ag n (2.13) 

a c " 

The practical computation of (2.13) is done using 0 *x 0 * terrestrial 
mean gravity anomalies Ag T for the numerical integration and a set of 
potential coefficients up to degree M as follows: 

N, = 4 ^ Jf S,(cosV)AgTd<r + ^ (Qln + W 1 n )Ag* (2.14) 

®c n 

To derive expressions for the global RMS undulation error associated 
with (2.14), we will follow the concepts given in (Christodoulidis, 1976): 
The global RMS undulation error will have four error components, due 
to: 

1. Erroneous Ag T (mean values of 0*x0* blocksize). This error will 
be called "propagation error". 

2. Finite blocksize (i.e. 0* is not infinitely small) used for the 
Stokes' numerical integration. This error will be called "discretion 
error". 

3. Erroneous potential coefficients. This error will be called 
"commission error". 

4. Limited degree of expansion of the potential coefficients used. 
This error will be called "ommission error". 

Notice that there is an analogy between error sources 1 and 3 and 2 
and 4. Errors 1 and 2 are to be interpreted as the errors due to 

A _ A 

erroneous input data in (2.14) (Ag T , Ag s ) whereas errors 2 and 4 are to 
be interpreted as the errors due to the limited degree of expansion of 
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the terrestrial data (AgJ are taken only to degree n = n T ) and the 
potential coefficient data (Agg are taken only up to degree n = M). 

Using the rule of thumb for the maximum degree n T of the 
terrestrial error degree variances, we have 

n T = 18O*/0* (2.15) 

where 0* is the blocksize (in degrees) of the mean gravity anomalies 
used in (2.14). 


We then furthermore assume that 

AgJ = hi + 4 0 * n * n T (2.16) 

and M[(AgJ) 2 ] = C n n T < n (2.17) 

with C n taken from a model of anomaly error degree variances. 


A common model is the Tscherning-Rapp model (Tscherning and 
Rapp, 1974): 


C 


n 


— D cn+2 

(n-2) (n+B) 


3 ^ n < ® 


(2.18) 


with A = 425.28 mgal 2 ; B = 24, S = 0.999617, C 2 = 7.5S 4 mgal 2 and will be 
used for the numerical applications. The C n values refer to a sphere of 
radius equal to the mean earth radius. 


Similarly, for the potential coefficient set, we assume: 


Agn = + eg 


0 - n - M 


(2.19) 
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M((AgJ*) a ] * C n M < n (2.20) 

Combining (2.13) and (2.14) together with the assumptions (2.16), (2.17) 
and (2.19), (2.20) we obtain the error when computing the geoid 

undulation through (2.14): 

5N 1 = N, - N* = + «N 1>2 + <SN 1>3 + 5N 1>4 (2.21) 


with 

*N|,» 

= 4^ JJ Si(cosi,)l\ld* 

ff c 

(2.22) 


5N jf2 

= 4^v JJ s i(cos V) X Agjdtr 

4nr „=„ T+ , 

(2.23) 


5N i>3 

= 2 ? X (Qin + W, n ) C B 
A 7 n = o 

(2.24) 


«n 1}4 

= 2 ; J n+ [ Qin + W,n)Ag " 

(2.25) 


Using (2.5), the errors £N{ 5N 1>2 can be rewritten as: 

<5N, , i = JJ s,(cost&) | T cJd<r - JJ S, (cosY 1 ) Z T cSd<r = 

^ n=o T <T~<r c n ° 

= JJ S^cosf) |*«5d<r - JJ s7(cosV) n |^Sda, or 

(T <T 

«N|,i= 2? I T (X in - Q in ) e J (2.26) 

’ n=0 

Here X< n are the Fourier coefficients of Sj(cosV'). Similarly, 

«N|,i= 5§ £ (X 1n -Q in )AgJ (2.27) 

‘•y n=n T +i 

According to our discussion above we have the following errors: 

Error 6N, }1 is the error due to erroneous terrestrial gravity data 
(propagation error) 
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Error <5N 1)2 is the error due to discretized data in the numerical 
integration (discretion error) 

Error <5N 1>3 is the error due to erroneous potential coefficients 
(commission error) 

Error 5N, j4 is the error due to the limited degree of expansion of the 
potential coefficients (omission error) 

The expected global mean square errors of each error source 
described above can now be found, using (2.26), (2.27), (2.24), (2.25): 


«Nf,t = M(«Nf, 1 ) = (t 

R 1 2 n T 

J? Z (X ln - Qln) 2 ^n 

*f 9 n = o 

(2.28) 

5N ? j2 = M( 6 N? >2 ) = (t 

c)’ l (X ln - Q, n ) 2 c n 

-r J njn +1 

(2.29) 

fiN ? >3 = M( 6 N? >3 ) = (t 

I (Qln + W ln ) 2 <5C n 

(2.30) 

fiN ? >4 = M(fiN? >4 ) = (7 

Tz\ ? + ^in) 2 ^n 

-7 y n-M+l 

(2.31) 


The total global mean square error, assuming uncorrelated error sources 
is 

6N? = 6N? tl + 5N? >2 + 6N? )3 + 6N? >4 (2.32) 

and the total global root mean square error is from (2.32): 

5N, = («N?,i + 5N? >2 + 6N? >3 + 5N? >4 )* (2.33) 

We thus have derived general expressions for the geoid undulation 
computation and it’s expected global RMS error when a modification of 
Stokes* function is attempted. Each of the following four methods, or 
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even any arbitrary modification of Stokes’ function which follows the 
principal (2.2) can be easily formulated according to (2.14) and (2.28) - 
(2.31). The derivation has been given step-by-step because of its 
importance for the error analysis discussion, and for the understanding 
of the sections 2.2 and 2.3 that follow. 

2.2 The Conventional Methods 
2.2.1 Stokes’ Method 

In Stokes’ method (i=l) no modification of the Stokes’ function is 
attempted. 

In this case: 

W i ( cosy) = 0 (2.34) 

which means that all the Fourier coefficients W jn are zero: 

W ln = 0 0 < n < • (2.35) 

and from (2.2) if follows that 

S,(cosY0 = S(cos V) (2.36) 

The Fourier coefficients of S^cos'f) can be taken from the expansion of 
the Stokes’ function in spherical harmonics: 

S,(cos VO = l p n ( cos y) (2.37) 

n=2 n-1 

Expansion (2.37) can be rewritten as 

S x (cos^O = 1 X ln P n (cos^) (2.38) 

n = o ^ 


Comparing (2*37) and (2.38) the Fourier coefficients Xi n of S^cos^) are 
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given by (cf. Heiskanen and Moritz, 1967, p. 97). 

Xi„ = 


0 

2 

n-1 


0 * n * 2 
2 < n < ® 


(2.39) 


Thus, the undulation in Stokes’ method is computed according to (2.14) 
as: 

Ni = J| S(cosV') Ag T d<r + g” (2.40) 

with the components of the global mean square error from (2.28) - (2.31) 
and (2.35), (2.39) 


II 

IS 

$ 

I* F(X ln - Q ln ) 2 <r n 

n - 0 

(propagation error) 

(2.41) 

AN?, 2 = | 


l 2 ® f 2 ) 2 

E , IrT-T ~~ 

n=n T +i vn 1 } 

(discretion error) 

(2.42) 

AN ?, 3 = | 

$ 

\ 2 M 

l Q?n<5C n 

n-o 

(commission error) 

(2.43) 

an?, 4 = 

ii. 

1 2 « 

l ^ QfnC„ 

} n = M+i 

(ommission error) 

(2.44) 


The total global RMS error is 

6Ni = (AN? , i + AN? , 2 + AN ?, 3 + AN?, 4 )* (2.45) 

The truncation coefficients Q ln can be numerically computed by a 
subroutine written by Paul (1973) based on an accurate recursive 


formula. 
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2.2.2 Meissl’s Method 

To avoid the discontinuity of the kernel (2.9), Meissl (1971) 
introduced the following modification: 


W 2 (cosy) = | 
It follows that 

s 2 (cosy) = { 

and 

s a (cosy) = [ 


s 0 o * y * y 0 
o y 0 < y * n 


s(cosy) - s 0 o * y ^ y 0 
s(cosy) y 0 < y ^ rr 


o o < f ^ y 0 

S 2 (cosy) y 0 < y ^ 7 T 


(2.46) 


(2.47) 


(2.48) 


Notice that the kernel S 2 (cos y) is now a continuous function of the 
spherical distance y everywhere in 0 * y ^ ?r. 


We now proceed for the computation of the Fourier coefficients X 2ri , 

W a „, Q2n* 


Introducing the following recursive formulas for P n (y), y = cosy 


(n+l)P n+1 (y) + nP n _! (y) = (2n+l)yP n (y) 

n 


i 

(2.49) 

(l-y a )P;(y) = n(P B -i (y) - yP„(y)) 

n 


i 

(2.50) 

yPn(y) - Pn-i(y) = nP n (y) 

n 


i 

(2.51) 

(2n+l)P„(y) = P;+i(y) - P^-ity) 

n 


i 

(2.52) 


(the prime denotes differentiation with respect to the argument y) 
and denoting with y 0 = cosy 0 , we take: 

w 2 n = f W 2 (y)P n (y)dy = f S 0 P n (y)dy = S 0 f P n (y)dy, or 

J - 1 •’-l i 


E 


Also Q 2n = P S 2 (y)P n (y)dy = f a °S 2 (y)P n (y)dy 

J _ 1 J —1 

= J a °S(y)P n (y)dy - S 0 J a °P n (y)dy 
From (2.54) and (2.52) we obtain 


Q2n ®in 


(P„+i(yo) - Pn-i(yo)) n * 1 


s 0 (yo+l) 


n = 0 


(2.54) 


(2.55) 


An alternative formula of Q 2n can be derived, using (2.55) and (2.49). 


Now P n +i(y) 


- nr„_. 


(2.56) 


Substituting (2.56) into (2.55) we have: 


®2n “ ®in + 


(Pn-i(yo) - y<>Pn(yo)) n * 1 


(2.57) 


s 0 (yo+i) 


n = 0 


Note that for n = 0 the Q 20 coefficient usually given in the literature 


(Rapp, 1985a) as 


Q 20 - S 0 (l y 0 ) 


(2.58) 


corresponds to our Q 20 + W 20 = -S 0 (y 0 +1) + 2S 0 = S 0 (l-yo) 
Finally, the X 2n coefficients are 


X 2n = {_ i S 2 (y)P n (y)dy = J_^ S(y)P n (y)dy - S 0 J_P n (y)dy or 


-2S 

X 2 „ = 0 

2 

n-1 
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Now we are ready to give the final formulas for Meissl’s method. The 
undulation is computed as 

N a = JJ S a (cos^)4g T d<r + (2.60) 

<r c 

with the components of the global mean square error from (2.28) - (2.31) 
and (2.53), (2.55), (2.59): 


m,i = | 

LSI 

l2yJ 

I (^2n Q2n) 2<y n 

n-o 

1 2 « f 9 12 

(propagation error) 

(2.61) 

6Nl,i = | 

y 

E + [£i - c n 

n = n T +i KXl 1 } 

(discretion error) 

(2.62) 

m,3 = 

( 2 ^) 2 [ n | a 0i" sc " + <«»• + 

W 20 ) 2 £C 0 ] 

(commission error) 

(2.63) 

5N|,4 = 1 

l2 r - 

\ l A QlnC n 

n-M+i 

(ommission error) 

(2.64) 


The total global IMS error is 

5N 2 = ( <5N! f i + 5Ni, 2 + tfNl, 3 + SNi,.)* (2.65) 

Expressions (2.61), (2.62), (2.63), (2.64) and (2.65) are equivalent to the 
expressions for the global mean square undulation error of Meissl’s 
method given in (Engelis, et al. f 1985a) if we are to start the summations 
from two instead of zero. 

2.3 The New Methods 
2.3.1 Molodenskii’s Method 


We will now present the Molodenskii’s modification to the original 



19 


Stokes’ formula (2.1) (Jekeli, 1980, p. 22). The function of -\ 
that is to be subtracted from the Stokes’ function S(cos^) is 


n 9k+l 

W^cos^ 1 ) = I — s- W 3k P k (cosV-) 

If — o * 


Using (2.2) we take 

n 2k+l 

S 3 (cos^) = S(cosV) “ I — o" W 3k P k (cos^) 

k— O * 


The Fourier coefficients X 3n , W 3n , and Q 3n can be immediately 
Rewritting (2.67) as 

S 3 (cosY0 = z x 3n P n (cosf) 

n — o ^ 


with 


l 3n 


'W 3n 

_2_ _ w 
n-1 " 3 " 

2 

n-1 


0 * n * 1 
2 ^ n * ii 
n < n < ® 


we obtain the Fourier coefficients X 3n . 


Furthermore, from (2.66) it is obvious that W 3n (0 * n * 
Fourier coefficients of W 3 (cosy). Finally, since 


S3 (cosio = ( ® s(cos1() 
we take (y = cos'Y, y 0 = cos^ 0 ) 

Q 3n = f S 3 (y)P n (y)dy = f y °S 3 (y)P„(y)dy = 

0 “I J - 1 

= J 9 ° (S (y) - l W 3k P k (y))P n (y)dy, or 

— 1 k = o z 


', W 3 (cosf) 
( 2 . 66 ) 

(2.67) 
computed. 

( 2 . 68 ) 

(2.69) 

n) are the 

(2.70) 
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fl ?lH-l fu n 

°3n = Qin ~ J J P k(y)P„(y)dy (2.71) 

The definite integrals from -1 to y, of the product of two Legendre 
polynomials P r (y)-P n (y) are analytically given in (ibid, p. 10). 

If 

erT1 ~ /_ 1 ^ r (y)‘P n (y)dy r, n * 0 (2.72) 

Then 

e rn = (r-n) (r+n+1) tnP r (yo)P n -i (yo) - rP n (y 0 )P r -i (yo) 

+ yo(r-n)P n (y 0 )P r (y 0 )] r * n r, n > 0 (2.73) 

Also, e r> © = 2 ^[ [P r +i(yo)P r -i(yo)] r > 0 (2.74) 

e 0 ,o = 1 + y 0 (2.75) 

e r,r = 2r+T £( 2 r_ l) e r-i,r-i + yo(P?(yo) + P?-i(y 0 )) - 2P r (y 0 )P r _ 1 (y 0 )] 

r > 0 (2.76) 

Formulas (2.73) - (2.76) are very useful for the numerical computation of 
e rn coefficients and can very easily be programmed. Notice that from 
(2.72) it holds 

e rn = e nr (2.77) 

Inserting (2.72) into (2.71) we have the truncation coefficients Q 3n : 
n 2k+l 

®3n — ^ln — i 2~ ^ 3 k e kn (2.78) 

In Molodenskii’s method the function W 3 (y) is to be defined as the 
"best" approximation to S(y) in the interval [-1, y 0 ], by minimizing the 
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norm of the difference: 

f 9 °[S(y) - W 3 (y)] a dy — » minimum (2.79) 

J — i 

After condition (2.79), the W 3n coefficients are computed (without the 
derivations, which can be found in (ibid, pp. 29-31)), from 

n 2r-+l 

W 3n = Z — s- u r h rn 0 < n * n (2.80) 

r=n *• 

with 

u r = f I ^ h rnQin 0 * r * n (2.81) 

k n=o ^ 

h rn = (2n+l)k-^| o (?)(?) (l-k) ,+l 0 * n < r 

h rr = 2r+l k_r r * 0 (2.82) 

p = r-n-1, q = r+n, k = cos 2 -^ 

It also follows immediately from condition (2.79) that 

Q 3n = 0 0 * n * n (2.83) 

The number n is the maximum number of "proper" harmonics that are 
removed from the Stokes’ function (see eq. (2.67)). For reasons that 
will be explained in Chapter 3, we choose 

n < M (2.84) 

Then, in Molodenskii’s method the undulation is computed as (ibid, eq. 
(95)): 

*3 = 4 ^- JJ S 3 (cosY0AgTd<7 + 2 ? Z W 3n Ag® + §- I Q 3 „Ag^ (2-85) 

' a ^7 n=o ^7 n-n+l 

w n 
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The components of the global mean square error from (2.28) - (2.31) and 
(2.69), (2.71), (2.80), and (2.83) are: 

< X >" - + - 0»n) a «.] 

(propagation error) 

5N| a = (o^l E (riT - Qsn] c n (discretion error) (2.87) 

'■*7' n=n T +i vn 1 ’ 

5N! 3 = (^) [ Z W§ n tfC„ + Z_ Q| n fiC n ] (commission error) (2.88) 

\£.y> i n =0 n=n+i J 

5N§ 4 = (^] Z Q§„C n (ommission error) (2.89) 

'■‘•I' n=H+l 

The total global RMS error is 

« i 3 = (fiNl.i + 5N|,2 + 6N?,, + «N§, 4 )* (2.90) 


2.3.2 Modified Sioberg's Method 

Proceeding as in Molodenskii’s method, we have that 
n 2k+l 

W 4 (cos D = z W 4k P k (cosT&) (2.91) 

k=o 

is the function to be subtracted from the Stokes’ function S(cos^). 
Then 

n 9k+1 

S 4 (cosV-) = S(cosf) - z W 4k P k (cosf) (2.92) 

k-o ^ 

or S 4 (cosf) = Z X 4n P n (cosV-) (2.93) 

n=o 

with 


X 


4n 


-W 4n 0<n<l 

- W 4 „ 2 < n * n 

2 

n-1 


(2.94) 


n < n < ® 
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Finally 


o o * 

S 4 (cOSY 0 ^ o < f i 7T 


Q«n = f S 4 (cosV)P n (cosf)d^ = f S 4 (y)p n (y)dy = 

J -7T J “I 

= f M °S 4 (y)P n (y)dy = f a °(S(y) - l W 4 k P k (y) )P n (y)dy or 

^ — l J — 1 k = 0 * 


n 2k+l 

Q*n ~ Qln - I o ^ 4 k e kn 

k — O L 


(2.96) 


with e kn from (2.73) - (2.76), 


Now, in general, Q 4n * 0, 0 * n < n, but W 4n = 0 for n < n * M 
since we again choose (see Chapter 3) n < M. Using the general 
formulation discussed in Section 2.1, we compute the undulation 
according to (eq. (2.14)) 

N* = A JJ S 4 (cos^)Ag T d<r + g £ (Q 4 „ + W 4n )Ag|* + -5 £ Q 4n Ag* 


2y n-n+i 


(2.97) 


with the components of the global mean square error from (2.28) - 


(2.31): 


*«.■ = (^nx <x *" - »♦.)*•. + 


(propagation error) (2.98) 


5N *,a " ( 27 ) n | n + ,(i^T “ Q «n) C n (discretion error) (2.99) 
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«N2,3 = (of) 2 [ I (Q4„ + w 4n )Mc n + ? Q2„«cJ 

’ v *7 / l n=o n=n+l J 

(commission error) ( 2 . 100 ) 

6 NJ 4 = [ 75 ^] Z Q 4 n C n (ommision error) (2.101) 

’ n=M+i 

and the total mean square error: 


m = 6Nl tl + 6Nl f2 + «NJ j3 + «N* )4 


( 2 . 102 ) 


In Meissl’s method care is taken for the kernel S 2 (cosf) to be 
continuous at i> - y 0 » in Molodenskii’s method the norm of the 
approximation of S(cos ip) through W 3 (cos^) is minimized (eq. (2.79)). In 
this method, the brilliant idea, due to Sjoberg (1986a) to directly 
minimize the expression of the mean square error ( 2 . 102 ) is utilized: 


— > minimum (2.103) 

To minimize (2.102) the partials are to be taken with respect to (up till 
now) arbitrary coefficients W 4n : 


= 0 k = 0, 1, ..., n 

dW 4k 


(2.104) 


Notice that the coefficients Q 4n are also linear functions of W 4n through 
eq. (2.96). The algebra will be given in Appendix A. We give the final 
result of the equations (2.104). They are equivalent to the linear 
system: 


Z a if 4 ^ h if * k 0, 1, ..., n 

t=o 


(2.105) 
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where a kt = a tk = (<r t + <5C k )<5 kt =- e kt (<r t + <5C t ) 


e kt (<r k + 6C k ) 


2k+l 2t+l f 8 


2 2 


[ E e tn e kn(®n ^n) 


+ E ®tn®kn( ff n ®n) "*■ JI e tn®kn^n "*■ ^n)] 
n=H+i n=n T +i 1 

9W+1 r M 

and h k = — s- I l e kn [Q ln (<r„ + <5C n ) - X ln <x n ] 

* l n = 0 


(2.106) 


+ D e knt®ln( ff n + ®n) Xi n ff n ] 
n = M+i 


E ®knf®in(^n + C n ) Xi n ff n ]l + Xj k ff k Qi k (<T k + <SC k ) 

n-n T +l 


(2.107) 


Kroenecker’s <5 k t is defined as <5 kt = 


The computation of W 4n coefficients is achieved by inverting 
the (symmetric) matrix A = i aick -li] and multiplying the result by 


the vector H = [hj ... h k ... h„] T : 
[W 41 W 42 ... W 4 "] t = A-*H 


(2.108) 


Once the W 4n coefficients are computed, the X 4n coefficients are 
computed backwards through (2.94), the Q 4n through (2.96) and finally 
the undulation through (2.97) and it’s global mean square error 
components through (2.98) - (2.102). 


CHAPTER III 


ERROR ANALYSIS 

In this chapter the numerical application of the formulas giving the 
global RMS undulation error for all the four methods will be examined: 
Formulas (2.41) - (2.44) for Stokes* method; (2.61) - (2.64) for MeissFs 
method; (2.85) - (2.89) for Molodenskii’s method; (2.98) - (2.102) for 
modified Sjoberg's method. Irrespectively of which method will be the 
most precise (i.e. will give the smallest global RMS undulation error), all 
the methods will be considered for the actual geoid computations. 

It is clear by examining the global RMS error formulas for all the 
methods discussed in the previous chapter that we need models for the 
error anomaly degree variances due to the errors of both the 0 # x0° 
terrestrial gravity anomalies (<r n ) and the potential coefficients (<5C n ) to 
numerically apply the formulas. Also, the anomaly degree variances (C n ) 
are required, and they are modeled as discussed before (eq. (2.18)) 
according to the Tscherning-Rapp model. 

The anomaly degree variances are used to compute the ommission 
error for the degrees M<n<® and the discretion error for degree n T <n<®. 
We here list some values of the anomaly degree variances C n of the 
Tscherning-Rapp model (Table 1), for high degrees. 
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Table 1. Anomaly Degree Variances of the Tscherning-Rapp Model (mgal 2 ). 


n 

C„(mgal 2 ) 

180 

1.955 

200 

1.766 

300 

1.173 

400 

0.862 

500 

0.671 

2000 

0.098 

5000 

0.012 


Tests using different values of A are later presented in this chapter. 


The error anomaly degree variances due to erroneous potential 
coefficients are taken from the errors 6C nm , 6S nm of the fully normalized 
potential coefficients of the OSU86D solution (Rapp and Cruz, 1986a). 
The above set of potential coefficients does not contain any zero or first 
order coefficients, and thus the summation when computing the 
commission error, will start from 2 

<5C n = r 2 (n-l) 2 I (5C nm 2 + 5S nm 2 ) 2*n*M (3.1) 

m— o 

Notice that in writing (3.1) we assume that the errors in the potential 
coefficients are uncorrelated. 


The numerical values of <5C n for some degrees are listed in Table 2. 
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Table 2. Error Anomaly Degree Variances due to the OSU86D Potential 
Coefficient Errors (mgal 2 ). 


n 

6C n (mgal 2 ) 

2 

0.644- 10 -6 

3 

0.745- 10" 4 

4 

0.296- 10“ 3 

5 

0.256- 10“ 2 

10 

0.568- lO -1 

20 

0.274 

100 

0.522 

150 

0.936 

180 

1.300 


For the error anomaly degree variances due to the errors in the 0*x0* 
terrestrial anomalies, various error models were examined. For all the 
models we assumed Agj = Ag J = 0, and thus for the computation of the 
propagation error (eq. (2.28)) the summation will again start from 2. 


MODEL A 

The error anomaly degree variances are taken from the numerical 
integration of the error covariance function C s (f) of the 6'xlO’ free-air 
anomalies in Europe (Weber and Wenzel, 1982). 


<r n = 


2n+l 

P n 2 


C e (Y0h(lMm a x )P n (cos^)sin^d^ 2*n*n T 


(3.2) 


FO 


where: 


Pn 

c e(V) 
r 

T^wax 


Pellinen smoothing factors for a spherical cap equivalent 
to a 6’xl0’ block 

22e -4 V'* + 3 mgal 2 error covariance function 
spherical separation in degrees 
maximum distance of integration = 10* 
step of integration = 1' 
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h(Y’»^max) = ^ [l + cos y^ ] Hanning smoothing window, to smooth out 
^ 1 r max'' 

the side lobes produced by the limited length (10*) of the 
error covariance function 

n T 1800 maximum degree at the anomaly error degree variances 

The above error covariance function implies an error variance C £ (0) = 25 
mgal 2 and a correlation length X = 0.*21. The correlation length X is 
defined as: 

C,(X) = *C £ (0) 

MODEL B 

Another model proposed by Sjoberg (1986b) is based on the 
following form of the error covariance function: 

°*<+) = c<1 - u) [ (l-2ucos^)» " 1 - ucos *] (3 - 3) 

where the coefficients u and c are computed by specifying the variance 
C £ (0) and the correlation length X. Then the error degree variances are 
given by the closed formula: 

<r n = c(l-u)u n 2^n<® (3.4) 

For this form three separate sets of C £ (0) and X were tested: 

MODEL B.l C £ (0) = 25 mgal 2 , X = 0.*21 (to comply with MODEL A) 

MODEL B.2 C £ (0) = 25 mgal 2 , X = 0.*1 

MODEL B.3 C £ (0) = 25 mgal 2 , X = l.*0 

The latter two models (B.2 and B.3) check the sensitivity of the global 
RMS undulation error with respect to the correlation length of the error 


covariance function. 
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MODEL C 

In case that we assume uncorrelated errors of "a” mgal between the 
0*x0* anomalies) then the error anomaly degree variances can be given 
by (Wichiencharoen, 1984): 

<r n = (2n+l) 0 2 2*n*n T (3.5) 

where 0 is the blocksize in radians; 

nx = 180/0* 0* being the blocksize in degrees. For this model we 

used 0* = 2/60 (2*x2' anomalies). 

Notice that the error degree variances given in (3.5) increase linearly 
with the degree. This does not happen with any of the models A, B.l, 
B.2, B.3. Numerical values of the error anomaly degree variances implied 
by models A, B.l, B.2, B.3 and C are listed in Table 3. 


Table 3. Error Anomaly Degree Variances (in mgal 2 ) Implied by Models A, 
B.l, B.2, B.3, C. 



MODEL ! 

n 

A 

B.l 

B.2 

B.3 

C 

2 

0.034 

0.053 

0.025 

0.257 

3. 38*10-® 

3 

0.048 


0.025 

0.255 

4.73*10“® 

4 

0.061 


0.025 

0.252 

6.05*10“® 

5 

0.073 

1 

0.025 

0.250 

7.40*10“® 

10 

0.121 

0.052 

0.025 

0.237 

1.42* 10“ s 

20 

0.133 

0.051 

0.025 

0.214 

2.75* 10“ s 

100 

0.034 

0.043 

0.023 

0.093 

1.35* 10“ 4 

1000 

0.010 

0.006 

0.009 

0.837* 10“ 5 

1.35* 10“ 3 

1500 

0.015 

0.002 

0.006 

0.472* 10“ 7 

2.02* 10“ 3 

1800 

0.037 

0.001 

0.004 

0.211 lO" 8 

2.43* 10“ 3 

5000 


0.128-10-* 

0.162*10“ 3 

0.861 • 10 -23 

6. 73 • 10“ 3 


For comparison purposes for all the four methods we will initially 
use Model A. Then, models B.l, B.2, B.3, C will be tested using the 
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modified Sjoberg's method only. 

The choice of n (see Chapter 2) to be less than the maximum degree 
of the potential coefficients used is discussed by Jekeli (1980, p. 33, fig. 
10). Using the Molodenskii’s method Jekeli showed that one obtains the 
lowest global RMS undulation error when n < M. However the above 
conclusion is based on errorless gravity data: 

<r n = 0 2*nO (3.6) 

Sjoberg (1986a, p. 14, fig. 2) also shows that the lowest global RMS 
undulation error is achieved if a "hybrid" solution is used with n = 20, 
M = 180. Finally, the choice n = M, in the case of the modified 

Sjoberg’s method with M = 180, will require the inversion of a symmetric 
matrix (see eq. (2.108)) of dimension 180x180 and the computation effort 
does not pay back in terms of gain in precision, as it will be shown 
later in this chapter. 

Figure 1 shows the global RMS undulation error for all the four 

methods for capsizes 0*, 1*, ... , 10*. The cases n = M = 20 are also 
shown for the Molodenskii’s and modified Sjoberg’s methods, which for 
i>o < 4* have large (greater than 70 cm) error due to the large 

ommission error. (For the rest of this discussion, when we use the 
terms "Molodenskii’s method", or "modified Sjoberg’s method" we will 
refer to the cases n = 20, M = 180, unless otherwise stated.) 

The gain in precision for all the capsizes when the modified 

Sjoberg’s method is used is clear. This results from the fact that the 




Figure 1. Global RMS Error in Undulation (in cm) for the Four Methods 
Using Model A for the Terrestrial Anomaly Errors and the 
OSU86D Potential Coefficient Errors. 


CAPSIZE (DEGREES) 
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modified Sjoberg’s method minimizes the global mean square undulation 
error itself (condition (2.103)). However, this improvement in precision 
for the capsize of 2* is on the order of 5 cm and may or may not be 
realized when actual geoid computations will be carried out using the 
above capsize. For capsizes the Molodenskii’s method is better 

than the Meissl’s method. For capsizes 3**^ 0 *6* the reverse happens. 
Then, for capsizes the Molodenskii’s method takes again the 

advantage. Stokes’ method compares well with the other methods for 
0*^ o < 2*. For larger capsizes the global RMS error increases faster for 
Stokes’ method than for any other method. This is due to the fact that 
the propagation error increases faster for Stokes’ method than for any 
other method as the capsize increases. Notice also from the above 
diagram that increasing the capsize does not substantially decrease the 
error in the case of the modified Sjoberg’s method, which is welcome 
when actual geoid computations are performed. 

In case that we assume <r„ = <r x - <5C 0 = 0 (see also Chapter 10), the 
modified Sjoberg’s method results in a singular matrix A to be inverted 
for =0*. This follows immediately from the orthogonality of the 
Legendre’s polynomials: (y 0 = cosO = 1) 

e rn = { P r (y)Pn(y)dy =0 r * n (3.7) 

' — 1 

Then the elements of matrix A from eq. (2.106) will be: 

a tr = 0 for t * r; a 00 = a cl = a 10 = 0 
and matrix A has two rows equal to zero. The value of the global RMS 
error for this case can be computed using the values of the commission 
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and omission error of the OSU86D set. (The propagation and the 
discretion error are zero). Also, matrix A is near-singular for the other 
nonzero capsizes (see also Colombo’s method discussed in (Jekeli, 1980, 
p. 38), where the condition fiN c = 0 results in an ill-conditioned system 
of equations). Further investigation is required that will result in some 
type of regularization of matrix A. In Sjoberg’s method the minimization 
of <5NJ is equivalent to the minimization of the weighted sum of the 
residuals V T PV = min in a traditional least squares adjustment. See also 
(Wenzel, 1982). 

The individual error components of the modified Sjoberg’s method 
are shown in Figure 2, and it can be seen that the major contibution to 
the error budget comes from the propagation and the commission error 
which for = 3* are on the order of 45 cm and 40 cm respectively. 
The omission and the discretion error are very small (less than 8 cm for 
i > o * 2*). 

We then tested models A-C for the error anomaly degree variances 
using the modified Sjoberg’s method. For these tests and for the rest 
of this discussion the capsize of 2* has been used, unless otherwise 
stated. The global RMS undulation errors due to the above models are 


shown in Table 4: 
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Figure 2. The Various Error Components of the Global RMS Undulation 
Error (in cm) for the Modified Sjoberg's Method (n □ 20, M = 
180). Model A for the Terrestrial Anomaly Errors and the 
OSU86D Potential Coefficient Errors Have Been Used. 
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Table 4. Global RMS Undulation Error (in cm) for the Modified Sjoberg’s 
Method Using Different <r„ -models for Anomaly Errors. 


MODEL 

6N A (cm) 

A 

60 

B.l 

58 

B.2 

53 

B. 3 

72 

C 

35 


We see by comparing model A and B.l that for the same error 
variance C e (0) and correlation length X, the form of the error covariance 
function does not give any significant difference for the total error 
budget. But the correlation length plays an important role and it can 
increase the error from 52.70 cm (X = O.'l) to 72.10 cm (X = 1*). 
Finally, assuming uncorrelated errors, the total error in undulation is 
reduced to 35.27 cm which is the most optimistic case. 

The choice of the error variance C c (0) is also a very important 
factor for the modified Sjoberg’s method. To demonstrate this we will 
show results of actual geoid computations in Maui Island (located in the 
vicinity of Hawaii). Using eq. (2.97) with n = 20, M = 180 we computed 
the geoid undulation for station 7210 in Maui from 2'x2' terrestrial and 
altimetric anomalies and the OSU86F set of potential coefficients. For 
more details see the next chapters. In Table 5 with N c we denote the 
first term in (2.97) (cap contribution) and with N p we denote the last 
two terms in (2.97) (i.e. contribution from the outer zone). Then N = N c 
+ N p . 
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Table 5. Actual Geoid Computations in Maui (Station 7210) with the 
Modified Sjoberg’s Method Using Different Error Variances. 



c»(0) (mgal 2 ) 1 


25 

10 

625 

N c (m) 

17.80 

19.17 

7.38 

N p (m) 

0.86 

-0.48 

9.43 

■m 

18.66 

18.69 

16.81 


Results of Table 5 are to be compared with the value of 20.49 m for 
the undulation of the station (See chapter 4). Notice that the smaller 
the error variance, the larger weight will be given to the terrestrial 
gravity data. For the pessimistic selection of the C e (0) to be 625 mgal 2 , 
the value of the undulation shows a discrepancy of ~3.50 m with the 
value of 20.49 m, whereas for the other two selections of C, ; (0) the 
discrepancy is ~1.80 m. 


Generally, each of the four methods we examine assigns different 
weights to the terrestrial gravity data and to the potential coefficients. 
To see this, we can rewrite (2.14) as 

Ni = 4^b \\\ Si (cosYOAg T dff - Jj Si (cosf)Ag T d<r] + ^ I (Q jn +W in )A| n 
7 <r <r-<r c 7 n 0 

=5? [ S T (X in -Qin)Ag l + S (Qin+W 1n )A| n 1 

^7 L n=o n -o J 

= 2 ? [ ^PiSAgJ + I P? n A| n ] (3.8) 

^7 L n~ O n — O J 

Thus, the weights that are used for the terrestrial data are: 


p. T = 
r i n 


p. T 
r l n 


p . T+p . s 
r i n r i n 


. £ i n 
*in+W in 


0^n-n T 


(3.9) 


For the potential coefficients they are: 
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P 


In 


PiS+P«8 



Qfn+Wfti 

X,„+W in 


O^n^M 


(3.10) 


Equations (3.9) and (3.10), in the case of Stokes' method (i=l),. are 
valid for n*2 and for Meissl's method (i=2) for n*l, so that the 

denominators of (3.9) and (3.10) do not become zero. The weights P^, 
Pjg for Stokes' method and P 2 J, P 2 g for Meissl's method are determined 
independently of the errors of Aj and Ag. (They are only functions of 
the capsize y 0 )* For the Molodehskii’s method P 3 J and P 3 g are 

determined immediately after the numbers n and M are given, but, again, 
without any refinement that would account for errors in A^, Ag. Only 
the modified Sjoberg's method gives weights that take into account the 
errors in Aj (through <t n ) and in Ag (through fiC n ). This is clearly 
shown in Table 5. The selection of the proper model for the anomaly 
degree variances however does not play an important role for the 

modified Sjoberg’s method, since this model is used to model the 

behavior of the anomaly spherical harmonics of high degrees only (n > 
M for the omission error and n > n* for the discretion error). To see 
this we list in Table 6 the various global RMS undulation errors for the 
modified Sjoberg's method, when different values of A for the C n model 
(see eq. (2.18)) are used: 


Table 6. Influence of the A-value of the C„ Model on the Global RMS 
Error for the Modified Sjoberg's Method. 


A (mgal a ) 

6N* (cm) 

80.00 

59 

425.28 

60 

10000.00 

65 
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The effect of different S-values (see eq. (2.18)) is expected to be 
neglibible for <5N 4 , but separate tests have not been carried out to show 
this. 


We now examine the global RMS undulation error for the 
Molodenskii’s and modified Sjoberg’s method in case that a higher 
degree of expansion M is used for the potential coefficients. Table 7 
clearly shows that increasing the degree of expansion from M = 180 to M 
= 250 the global RMS error remains essentially the same. 

Table 7. Global RMS Undulation Error (in cm) for Molodenskii’s and 
Modified Sjoberg’s Methods Using Two Different Degrees of 
Expansions. 



<5N 

(cm) 

METHOD 

M = 180 

M = 250 

Molodenskii’s 

61 

61 

Modified Sjoberg’s 

60 

60 


Finally, different values of n (number of harmonics in W 3 (Y I ) and 
W 4 (^)) influence the global RMS error only in the case of Molodenskii’s 
method as Table 8 shows (M = 180). 


Table 8. Global RMS Undulation Error (in cm) for Molodenskii’s and 
Modified Sjoberg’s Method Using Various Numbers of Harmonics 
in W,(f). 


Tl 

5N 3 

<5N 4 

5 

63 

61 

20 

61 

60 

40 

62 

60 

180 

71 

60 
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Table 8 shows that the combination n = 20, M = 180 gives the 
smallest global RMS error in the undulation in the case of Molodenskii’s 
method. This result is in analogy with (Jekeli, 1980, p. 33, fig. 10). 
The improvement we obtain in the case of the modified Sjoberg’s method 
(1 cm) is negligible and the computation effort much larger when we 
increase n. 

Summarizing, we see that the selection of the cr n and 6C n models is 
critical for the modified Sjoberg’s method, whereas the C„ model does 
not play an important role. Although we can have a reasonable <5C n 
model using eq. (3.1), it is very difficult to have a reasonable <j n model 
for the terrestrial gravity anomalies, because for the majority of the 
cases, either 

a) No error estimates are available, or 

b) Unrealistic error estimates are given. 

Even if reasonable error estimates were available for the given gravity 
data, the algorithms to compute the error covariance function 
and/or the corresponding <r n model have not yet been developed 
satisfactorily: The method of computing an error covariance function by 

comparing two independent data bases for the same region (Weber and 
Wenzel, 1982) suffers from the fact that two independent data sets are 
rarely given for the same region since the gravity measurements are 
usually made once within a broad time period. On the other hand, the 
empirical error covariance function computation based on actual error 
estimates will not give a unique solution due to the plus or minus sign 


of the error estimates. 
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Finally, the selected <x n models are only approximate and the 
assumption that they are valid everywhere on the earth is not 
rigorously true. Also, the reference surface that is associated with 

these models has been taken to be the sphere of radius R, or, in a first 
approximation, the geoid. Changing the radius of the reference surface, 
by several hundred meters, will cause a small change in the cr n models 
and such a change will not significantly affect the error analysis. 
However from the theoretical point of view, it is critical to associate a 
surface with the selected models. 


CHAPTER IV 


LASER STATIONS 

The coordinate system that was selected as the consistent 
coordinate system for defining the laser station coordinates was the SL6 
system (Robbins et al., 1985). Although the laser station coordinates 
were also given in the SL5.1AP system (Smith et al., 1985), we selected 
the SL6 system as the most recent one. Initially, the geodetic 
coordinates (♦,X,h) of 46 laser stations referred to the SL6 solution were 
available from Robbins et al., (1985). The station number, name and 
occupation of these stations are given in Table 9 (Alfano, 1986). 
Stations 7907 in Peru and 7929 in Brasil were immediately excluded from 
the list of the stations that we would finally use for geoid undulation 
computations due to the lack of local terrestrial gravity information. We 
thus retained 44 stations whose coordinates had to be transformed from 
the SL6 reference ellipsoid (a e = 6378144.11 m; f = 1/298.255) to our 
adopted reference ellipsoid (we will call it OSU GRS: Ohio State 
University Geodetic Reference System) with parameters: 
a e = 6378136.0 m 
f = 1/298.257222101 

(4.1) 

GM = 398600.440 km 3 /sec 3 
a - 729115 -10" 11 rad/sec 

The three last parameters in (4.1) are the same as the parameters 
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NUMBER I NAME I OCCUPATION 


QUILAS Quincy, California 

QUILAS Quincy, California 

QUILAS Quincy, California 

SANDIE Otay Mountain, San Diego, California 

7110 MONLAS Mount Laguna, California 

7082 BEARLK Bear Lake, Utah 

7084 OVRLAS Owens Valley Radio Observatory, Big Pine, Calif. 

7114 OVRLAS Owens Valley Radio Observatory, Big Pine, Calif. 

7085 GOLDLS Goldstone, California 

7115 GOLLAS Goldstone, California 

7086 FTDAVS McDonald Observatory, Fort Davis, Texas 

7885 FORLAS McDonald Observatory, Fort Davis, Texas 

7112 PLALAS Platteville, Colorado 

7887 VANLAS Vandenberg Air Force Base, California 

7888 HOPLAS Mount Hopkins, Arizona 

7921 HOPLAS Mount Hopkins, Arizona 

7894 YUMLAS Yuma Proving Grounds, Arizona 

7063 STALAS GORF, GSFC, Greenbelt, Maryland 

7064 GSFCLS GORF, GSFC, Greenbelt, Maryland 

7100 GSF100 GORF, GSFC, Greenbelt, Maryland 

7101 GSF101 GORF, GSFC, Greenbelt, Maryland 

7102 GSF102 GORF, GSFC, Greenbelt, Maryland 

7103 GSF103 GORF, GSFC, Greenbelt, Maryland 

7104 GSF104 GORF, GSFC, Greenbelt, Maryland 

7105 GSF105 GORF, GSFC, Greenbelt, Maryland 

7069 RAMLAS Patrick Air Force Base, Florida 

7091 HAYLAS Haystack Observatory, Westford, Massachussets 

7120 MUILAS Lure Observatory, Mount Haleakala, Maui, Hawaii 

7210 MAULAS Lure Observatory, Mount Haleakala, Maui, Hawaii 

7090 YARLAS Yarragadee, Australia 

7943 ORRLAS Orroral Valley, Australia 

7067 BDILAS Bermuda Island 

7805 FINLAS Metsahovi, Finland 

7939 MATLAS Matera, Italy 

7835 GRALAS Grasse, France 

7840 RGOLAS Royal Greenwich Observatory, Great Britain 

8833 KOOLAS Kootwijk Observatory, Apeldoorn, Netherlands 

7834 WETLAS Wettzell, Federal Republic of Germany 

7838 SHOLAS Simosato Hydrographic Observatory, Japan 

7061 EASTER Easter Island, Chile 

7068 GRKLAS Grand Turk Island 

7092 KWJLAS Kwajalein, Marshall Islands 

7096 SAMLAS American Samoa 

7121 HUANIL Huahine, Society Island, French Polynesia 

7907 ARELAS Arequipa, Peru 

7929 NATLAS Natal, Brazil 
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of the GRS’80 reference ellipsoid. For the transformation of coordinates 
we followed the scheme (Rapp, 1985c): 

(♦j A,h) sl 6 — » (x,y,z) (4>,A,h) 0SU 

Provided the orthometric height H of the laser station is available 
the value of the undulation N T of the station is the difference 

N T = h - H (4.2) 

where h is the given ellipsoidal height. The reference point for the 
ellipsoidal heights of the laser station was given either as a fixed 
marker on the ground (Bench Mark) or as the intersection of the 
horizontal and vertical axes of the laser instrument. Accordingly, the 
reference point for the orthometric height should be the same, so that 
the value of the undulation computed from equation (4.2) is meaningful 
and useful for comparisons. In case that the reference point is not the 
same then discrepancies up to 3.5 m (maximum vertical separation 
between the Bench Mark and the instrument’s axis) can occur. Using 
the calibration information for the laser stations included in (Noll, 1983) 
we transformed the orthometric heights to refer either to the fixed 
marker or to the intersection of the instrument axis so that the 
orthometric and the ellipsoidal height would have the same reference 
point. This transformation was not possible for 11 stations because no 
calibration data (i.e. North-East-Up coordinates of the intersection of the 
instrument axis with respect to the fixed marker) was available for those 
stations. For the above 11 stations the orthometric heights were taken 
from (Alfano, 1986), and they are thought to refer to the intersection of 


the instrument axes. 
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In Table 10 we see the coordinates of the stations transformed to 
the OSU GRS, the orthometric height, the reference point (BM = fixed 
marker, IA = intersection of the laser instrument axes, UN = unknown) 
and the undulation computed using eq. (4.2) for the 44 laser stations. 
There are laser stations very close to each other (see for example the 
coordinates of the stations 7051, 7109 and 7886 in Quincy, California, in 
Table 10). It is useful to group those stations under a common name. 
This is shown in Table 11, where the laser stations with distances less 
than 100 km are grouped under the same name. 


Table 11. Grouping the Laser Stations with Distances Less than 100 km 


NAME 

STATIONS INCLUDED 

QUI 

7051,7109,7886 

SAN 

7062,7110 

OVR 

7084,7114. 

GOL 

7085,7115 

FTD 

7086,7885 

HOP 

7888,7921 

STA 

7063,7064,7100,7101,7102,7103,7104,7105 

HAW 

7120,7210 


The distribution of the 44 laser stations is shown in Figure 3 and 
the distribution of the 11 laser stations with unknown height references 
is shown in Figure 4. The accuracy of the undulation computed from 
equation (4.2) has not been rigorously computed, since no accuracies 
were available for the orthometric heights nor for the ellipsoidal 


heights. 
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Table 10. 


Laser Station Coordinates Referred 
6378136.0 in, f = 1/298.257222101) 


to the OSU GRS (a e = 


NUM 

NAME 

DD 

LATITUDE 
MM SS.SSSS 

LONGITUDE 
DD MM SS.SSSS 

ELL. (M) 

ORT. (M) 

REF 

UND(M) 

705 1 

QU I LAS 

39 

58 

24.5651 

239 

3 

37.5530 

1060.9132 

1083 . 4800 

BM 

-22.5668 

7103 

QUILAS 

39 

58 

29.9961 

239 

3 

13.9490 

1 107.2592 

1129.8520 

BM 

-22.5928 

78 36 

GUI LAS 

39 

58 

30.0121 

239 

3 

18.0180 

1 1 10.5048 

1 129.8000 

UN 

-19.2952 

7062 

SAND IE 

32 

36 

2.6526 

243 

9 

32.7810 

989.5260 

1022.0000 

BM 

-32.4740 

71 10 

MONLAS 

32 

53 

29.9965 

243 

34 

38.2580 

1839.9153 

1870.7780 

8M 

-30.8627 

7032 

BEARIK 

41 

56 

0.8900 

248 

34 

45.5370 

1963.9328 

1976.5150 

BM 

-12.5822 

703 * 

OVRLAS 

37 

13 

55.6502 

241 

42 

15 . 1130 

1 179.0608 

1203.9580 

BM 

-24.8972 

7114 

OVRLAS 

37 

13 

57.2062 

241 

42 

22.2150 

1178.9648 

1203.8060 

BM 

-24.8412 

7035 

GCLPLS 

35 

25 

27.9573 

243 

6 

48.9170 

966.3705 

996.0500 

BM 

-29.6795 

7115 

GOLLAS 

35 

1 4 

53.8944 

243 

12 

28.9490 

1039.5651 

1069. 1450 

BM 

-29.5799 

7036 

FTDAVS 

30 

40 

37.2987 

255 

59 

2.4810 

1962.3776 

1983.1610 

BM 

-20.7834 

7335 

FORLAS 

30 

4 0 

37 . 3007 

255 

59 

2 . 4780 

1962.3310 

1983. 1610 

BM 

-20.8300 

7112 

PLALAS 

40 

10 

57.995 1 

255 

16 

26 .3360 

1502.5151 

1519.9130 

BM 

-17.3979 

7387 

VANIAS 

34 

33 

53.3514 

239 

29 

57.9780 

605.2623 

636 . 4000 

UN 

-31 . 1377 

7338 

HOPLAS 

31 

41 

6 . 3096 

249 

7 

l 8 . 5000 

2335.6675 

2337.8630 

UN 

-2.2015 

7321 

HOPLAS 

31 

41 

3.2166 

249 

7 

18.8370 

2353.9135 

2383.3870 

UN 

-29 .4735 

7334 

YUMLAS 

3 2- 

56 

20.9235 

245 

47 

48 . 6070 

242.6696 

271 .3000 

UN 

-23.6304 

7063 

STALAS 

39 

1 

13.3561 

283 

10 

1 9.7950 

20.2032 

53.1200 

BM 

-32.9168 

7064 

GSFCLS 

39 

1 

15.0981 

283 

10 

18.6050 

18. 1892 

51 .8260 

BM 

-33.6368 

7 100 

GSF 100 

39 

1 

15.4451 

233 

10 

47.6350 

1 1 . 1462 

43.6440 

BM 

-32.4978 

7101 

GSF 101 

39 

1 

16.1991 

283 

10 

42.8350 

9.3502 

42 . 4330 

BM 

-33.0828 

7102 

GSF 102 

39 

1 

1 4.374 1 

233 

10 

18.7920 

18.9272 

51.8150 

BM 

-32.8878 

7103 

GSF 103 

39 

1 

14.601 1 

283 

10 

18.7950 

18.8632 

51 .8030 

BM 

-32.9338 

7104 

GSF 104 

39 

1 

17.0761 

283 

10 

36.8380 

10.9342 

43.8760 

BM 

-32.9418 

7105 

GSF 105 

39 

1 

14.1531 

283 

10 

20. 1580 

20. 1202 

53.0530 

BM 

-32.9328 

7069 

PAMLAS 

28 

13 

40.6470 

279 

23 

39.2980 

-22.6006 

6.5290 

BM 

-29. 1236 

7091 

HAYLAS 

42 

37 

2 1 .6830 

288 

30 

44 . 3390 

92.9496 

120.5020 

BM 

-27.5524 

7 120 

MU I LAS 

20 

42 

27.3830 

203 

44 

38. 1020 

3068.7162 

3048.2530 

BM 

20.4632 

7210 

MAULAS 

20 

42 

25.9920 

203 

44 

3 8.6000 

3069.2871 

3048.7980 

IA 

20.4331 

7090 

YARLAS 

-29 

-2 

-47.4048 

1 15 

20 

43 . 107 1 

242.2921 

266.5630 

BM 

-24.2759 

7943 

ORR LAS 

-35 

-37 

-29.7502 

143 

57 

17 . 1240 

949.8848 

929.5000 

UN 

20.3848 

7067 

BD I LAS 

32 

21 

13.7566 

295 

20 

37.9270 

-22.0603 

9.7510 

BM 

-31 .8113 

7805 

F INLAS 

60 

13 

2.2828 

24 

23 

40.21 10 

79. 1806 

59.2300 

IA 

19.9506 

7939 

MATLAS 

40 

38 

55.7871 

16 

42 

16.6860 

536.9066 

490.52U0 

UN 

46.3966 

7335 

GRALAS 

43 

45 

16.8730 

6 

55 

1 5.8640 

1323.9484 

1281 . 1000 

UN 

42.8484 

7 3 40 

FGOLAS 

50 

52 

2.5551 

0 

20 

9 . 8620 

76.2630 

30.6300 

UN 

45.5830 

8333 

KOOLAS 

52 

10 

42.2392 

5 

48 

35. 1190 

94.4558 

49.8000 

UN 

44.6558 

7834 

VETLAS 

49 

8 

41 .7711 

12 

52 

40.9670 

662.0941 

614 . 4400 

IA 

47.6541 

7338 

SHOLAS 

33 

34 

39.7155 

135 

56 

13. 1890 

102.3637 

62 . 4000 

UN 

39.9697 

7061 

EASTER 

-27 

-8 
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Figure 3. Distribution of the 44 Laser Stations, With Known (33 Stations) and Unknown (11 
Stations) Height References. 
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Figure 4. Distribution of the 11 Laser Stations with Unknown Height References. 





CHAPTER V 


GRAVITY DATA PROCESSING 


5.1 Pre-processing of the Data 

The terrestrial gravity sources available by region are listed in 
Table 12 below (Ag = free-air gravity anomaly, Ag = mean Ag). The 
third column in Table 12 shows the number of laser stations in the 
corresponding region. 

Table 12. Gravity Sources Available by Region. 


REGION 

SOURCE 

# of LASER STATIONS 

UNITED STATES 

Point Ag 

27 

BERMUDA 

Point Ag 

1 

AUSTRALIA 

Point Ag 

2 . 

EUROPE 

6'xlO' Ag 

6 

JAPAN 

10'xlO' Ag 

1 


The gravity sources in Australia, Europe, and Japan are discussed in 
(Despotakis, 1986). The point gravity anomalies in the United States, 
Hawaii and Bermuda have been received in 1983 from the National 
Geodetic Survey; the 10'xlO' mean anomalies are also described in 
(Ganeko, 1982); and the 6'xlO' mean anomalies are also described in 
(Torge et al., 1983). These gravity sources cover only 39 laser stations 
as Table 12 indicates. For five (7061,7068,7092,7096,7121) oceanic stations 
the adjusted GEOS-3 and SEASAT oceanic altimeter data (Liang, 1983) 
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was used. For these stations an approximate geoid undulation was 
computed from the altimeter sea surface heights using a standard least 
squares collocation method using procedures similar to (Rapp, 1985b). 

Numerical details will be given in Chapter 8. 

Based on the error analysis results (Chapter 3, Fig. 1) the capsize 
for the cap contribution computations (first integral term in eq. (2.14)) 
was selected to be 2°. Thus, within a circle of radius 2° (with the 
center of the circle at the laser station) sufficient gravity data had to 
exist. More precisely, within the above circle we had to have full 
coverage of mean gravity anomalies so that the numerical integration of 
equation (2.14) would be possible. The gravity anomaly to be used in 
(2.14) has to: 

a) Refer to the reference formula of normal gravity implied by the 

constants of the adopted reference ellipsoid (see Chapter 4, eq. 

( 4 . 1 )); 

b) Be corrected for atmospheric effect; 

c) If possible, to be corrected for the terrain efffect. The 

undulation then has to be corrected for the indirect effect (see 
Section 5.2). 

All the given sources (Table 12) had their free-air anomalies 
referred to the GRS’67 reference formula. The conversion from this 
formula to the formula implied by the constants of the OSU GRS will be 
derived if we take into account that: 
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9<v 

A &1967 = Sobs ~ 71967 + h 

A 8osu = gobs " rosu + h 


(5.1) 

(5.2) 


where Ag, 96 7 : free-air gravity anomaly referred to the GRS’67 
Agosu : free-air gravity anomaly referred to the OSU GRS 
7 1 9 6 7 : normal gravity on the GRS’ 67 ellipsoid 
7o su : normal gravity on the OSU GRS ellipsoid 

3 y 

: free-air gradient of the normal gravity 
h : elevation of the station 


The quantities 9y/ 3h-h and g obs are the same for both systems. Thus, 
subtracting (5.1) from (5.2): 

Ag0SU = Agi 9 67 + 7,967 - 7osg (5.3) 

Knowing that (International Association of Geodesy, 1971) 

71967 = 978031.85(1 + a'sin 2 ^ + b'sin 2 2$) mgal (5.4) 

with a' = 0.0053024 and b' = -0.0000059, 

we can write (Heiskannen and Moritz, 1967, eqs. (2-116) and (2.105a)) 

7osu = (l - | m - jf e * 2 m) (1 + a'sin 2 * + b'sin 2 2+) (5.5) 

assuming that a' and b' will be the same as in (5.4) which is correct 
for the decimal digits we give. We also have: 


and a,f,GM,u from (4.1). 

We then obtain 

7 osu = 978032.83(1 + a'sin 2 $ + b'sin 2 2$) mgal (5.6) 

From (5.3), (5.4) and (5.6), we take: 

Agosu = Ag , 967 + (978031.85 - 978032.83) (1 + 0.0053024sin 2 * 

- 0.0000059sin 2 2$) mgal, or 


AgoSU ~ Ag, 967 


0.98 mgal 


(5.7) 
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Thus, from the anomalies referred to the GRS'67 reference formula, we 
have to subtract 0.98 mgal, to convert them to the OSU GRS reference 
formula. 

The next step is to apply the atmospheric correction to (5.7) 
considering the polynomial expansion of the atmospheric correction given 
in Wichiencharoen (1982a): 

fig A = 0.8658 - 9.727* 10 -5 *H + 3.482* 10“ 9 *H 2 (5.8) 

where H is the elevation of the station, and fig A , in mgal, is the 
atmospheric correction computed truncating (5.8) to the second term and 
giving an approximate value H for the elevation of the station. Then 
the free-air anomalies to be used in equation (2.14) were computed as: 

Ag = Agosu + (5.9) 

The above procedure applied to both point and mean Ag given from the 
various sources listed in Table 12. 

To determine the coverage of the gravity anomalies surrounding the 
laser stations, windows of A4 - 5* in latitude and AX = INTEGER(5*/cos$) 
+ 1* in longitude (♦: latitude of the laser station) centered at the laser 
stations were used. Full coverage and mean anomalies (as it was the 
case with the anomalies surrounding stations 7838 in Japan and 8833, 
7834 in Europe) did not require any further process. Insufficient 
coverage and mean anomalies (as it was the case with the anomalies 
surrounding stations 7805,7939,7835 and 7840 in Europe) required that 
we fill-in the gaps. For this purpose we used the OSU86F (Rapp and 
Cruz, 1986b) potential coefficient set up to degree 360 to compute the 
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missing Ag as the point Ag evaluated at the centers of the 6'xlO' empty 
cells through the formula: 

rjM 360 n n 

Ag = —2 l (n-l)l-l i (C nm cosmX + S nm sinmX)P nm (sinf) (5.10) 

r n= 2 vr ' m=o 

where: GM,a : the constants of the OSU GRS given in (4.1) 

C n m>S nm : the fully normalized potential coefficients of the OSU86F 

solution; the C 20 i C 40 and C 60 refer to the OSU GRS 
P nm (sin$): the fully normalized associated Legendre functions; $ is 

the geocentric latitude of the computation point 
r : the geocentric radius to the computation point 

X : the longitude of the computation point. 

An example of a case with insufficient coverage of the 6'xlO' Ag is 
shown in Figure 5 for station 7805 in Finland. The dots represent 
available 6'xlO' Ag and the dashes represent fill-in anomalies using the 
OSU86F field as described above. 

The point gravity anomalies surrounding the laser stations in the 
United States, Hawaii, Bermuda and Australia were transformed to 2' x2' 
mean anomalies using interpolation techniques (collocation or weighted 
average) that are discussed in Section 5.2. For the application of the 
collocation method using the five closest points it is necessary that the 
matrix (C 5 j + Djj) is non-singular (see Section 5.2, eq. (5.11)). If the 
five closest points are very close to each other, which was the case for 
most of the distributions of the given point anomalies, then the above 
matrix becomes singular. To avoid this, a thinning procedure was 
applied to the originally available point anomalies. This procedure was 


LONGITUDE 


Figure 5. Example of Insufficient Coverage of the 6'xlO' Ag 
Surrounding Station 7805 _in Finland. Dots Represent 
Originally Available 6'xlO' Ag and Dashes Represent Fill-in 
6'xlO' Ag Using the OSU86F Potential Coefficient Set Up to 
Degree 360. 
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as follows: For each 2'x2' cell in which the whole area was subdivided; 

one value, the one which was closest to the center of the 2'x2' cell was 
kept. The thinned data set created this way was then examined to 
ensure sufficient coverage. If the coverage was sufficient, the 
prediction of the 2'x2' Ag was done as discussed in Section 5.2. An 
example of such coverage for station 7082 in Bear Lake, Utah is shown 
on Figure 6. The dots represent the point anomalies of the thinned 
data set. If the coverage was not sufficient (i.e. areas larger than 
l’xl* were completely empty), then 

a) For coastal stations or stations on islands the 0.*125 x 0.*125 
altimeter anomalies of the world data base described in (Rapp, 
1985b) was used to fill-in the empty areas. Examples of this case 
are shown in Figure 7 for a coastal station (station 7887 in 
Vandenberg, California) and in Figure 8 for the stations 7120,7210 
on the island of Maui. The dots represent point Ag. 

b) For continental stations, for which no altimeter anomalies were 
available, the fill-in anomalies were generated at the centers of 
the missing 2' x2' cells using eq. (5.10) as described above. 
Figure 9 illustrates an example of such a case for stations 7888, 
and 7921 in Mount Hopkins, Arizona. The dots represent point 
Ag. 

After the fill-in procedures, the 2'x2' Ag were predicted as discussed in 
Section 5.2. The thinning procedure was not applied to the data 
surrounding station 7067 in Bermuda (see Section 5.3). 
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Figure 7. Example of Insufficient Coverage of Point Ag at the Coastline 
Station 7887 in Vandenberg, California. The Dots Represent 
Point Ag. 
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Figure 8 


Example of Insufficient Coverage of Point Ag at the Two 
Stations 7120, 7210 in the Vicinity of Hawaii. The Dots 
Represent Point Ag. 
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Some statistics concerning the terrestrial gravity data processing 
are shown in Table 13, The first column is the number or the name of 
the region (in case of groups of stations within the same region, see 
Chapter 4). The second column is the AX dimension selected for each 
window surrounding the laser station. The third column shows the 
number of point or mean Ag available for each region before ((1)) and 
after ((2)) the thinning. If no thinning was applied to the data (e.g. 
6'xlCT Ag data) this column is left blank (" — M ). The fourth column 
represents the fill-in procedure chosen: A M — " means that no fill-in 
anomalies were used; "P" means that the fill-in anomalies were taken 
from the OSU86P potential coefficient set and "A" means that the 0*.125 
x 0 # .125 altimeter anomalies were used to fill-in the gaps. The fifth 
column gives the final gridded form (or equivalently: the blocksize of 
the mean anomalies) in which the original data were transformed to for 
the numerical integration. The sixth and seventh columns give the 
variance and the arithmetic mean of the final mean anomalies over the 
sample area. 

For each grid cell the anomaly was computed from (5.10) using the 
OSU86F potential coefficients up to degree 180 and then the arithmetic 
mean over the sample area of those anomalies was computed. This mean 
is shown in column eight. Finally, column nine gives the mean elevation 
that was used for each particular window for the atmospheric effect 
computation through equation (5.8). 



61 


Table 13. Information Concerning the Terrestrial Gravity Data 
Processing. 


# OR 
NAME 




GRIDDED 

FORM 

VAR 

[mgal 2 ] 

X 

[mgal] 

[mgaS] _ 

H 

[m] 

mam 

tali 

QUI 

6 

43915 


— 

2'x2' 

1902.29 

7.93 

6.13 

500 

SAN 

6 

62506 


p 

2x2' 

1231.22 

-11.21 

- 9.71 

500 


7 

18654 

9882 

— 

2'x2' 

1901.87 

5.55 

16.26 

1000 

OVR 

6 

39966 


— 

2x2' 

2191.67 

2.84 

2.01 

700 

GOL 

6 

52428 


— 

2'x2' 

2288.04 

- 6.32 

- 5.70 

700 

FTD 

6 

12166 

5961 

— 

2'x2' 

518.89 

- 0.38 

4.58 

1000 

7112 

7 

26953 

11150 

— 

2'x2' 

1353.25 

16.75 

16.53 

1000 

7887 

7 

96744 

14345 

A 

2'x2' 

1462.18 

-13.6 

-14.15 

500 

HOP 

6 

13343 

5447 

P 

2'x2' 

435.85 

- 9.07 

- 7.08 

1000 

7894 

6 

23765 

7296 

P 

2'x2' 

849.12 

-10.94 

- 7.28 

100 

STA 

6 

29561 

11705 

— 

2x2' 

513.82 

- 3.21 

- 1.09 

300 

7069 

5 

43186 

5986 

A 

2'x2' 

415.67 

5.57 

2.79 

100 

7091 

7 

31095 

10357 

A 

2'x2' 

336.75 

0.89 

2.84 

100 

HAW 

5 

715 

458 

A 

2x2' 

10802.81 

25.30 

23.79 

100 

7067 

6 

15731 

— 

A 

2'x2' 

705.18 

-11.00 

-11.33 

0 

7805 

9 

1859 

— 

P 

6'xlO' 

293.95 

-15.95 

-15.71 1 

500 

7939 

7 

1886 

— 

P 

6'xlO' 

2927.86 

15.33 

15.39 

500 

7835 

6 

1747 

— 

P 

6'xlO' 

2429.60 

9.29 

8.89 

1000 

7840 

7 

1804 

— 

P 

6'xlO' 

128.03 

- 4.16 

- 1.87 

300 

8833 

8 

2400 

— 

— 

6'xlO' 

344.54 

- 0.68 

1.81 

100 

7834 

8 

2400 

— 

— 

6'xlO' 

686.79 

21.91 

21.32 

300 

7838 

6 

1080 

— 

— 

10x10' 

1550.13 

22.87 

25.42 

— 

7090 

6 

11835 

4091 

A 

2'x2' 

920.58 

-11.85 

-12.75 

100 

7943 

6 

9804 

4730 

A 

2'x2' 

1645.25 

14.01 

17.33 

500 
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Some blunders ( e.g . point anomalies of 1700 mgals) were found for 
some regions in the western United States, but they were easily 
removed by applying an editing method (i.e. delete the values of the 
point anomalies if they are outside certain specified limits). Thus the 
conclusion is that editing of the data is necessary, especially when we 
deal with large amounts of data, to avoid any blunders that the data 
may contain. 

5.2 Gravity Predictions 

For the transformation of the randomly distributed data points (see 
Figure 9 for example) into a gridded form of 2'x2' Ag, predictions of 
the 2'x2' Ag had to be carried out based on the thinned data covering 
the areas surrounding the laser stations and merged whenever 
necessary with either the altimeter data or the OSU86F anomalies, as 
discussed in the previous Section. Two basic types of prediction 
methods were considered (Cruz, 1983): 

a) Collocation method : The 2'x2' Ag were computed as the point Ag at 

the centers of the 2'x2' cells based on the five closest Ag to the Ag to 
be predicted as follows: 

Ag = + 0ij) -1 (Ag 1 - Ag) + Ag (5.11) 

where: A^ is the predicted point value at the center of the 2' x2' cell 

(=2'x2' Ag) 

Agj is the 5x1 vector of the closest five anomalies to the 
predicted point 

Ag = I Ag,/5 

i = i 
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Djj is a 5x5 diagonal matrix of the noise of the five data 
points 

C p j = Cpj(^): 1x5 covariance matrix between the prediction 
point and the five data points 

Ctj = C 1 j (V) '• 5x5 variance covariance matrix between the five 
data points 

y/n spherical separation. 

To compute matrices C p j (Y>) and C f j(V) the covariance function based on 
the Tscherning-Rapp model for the anomaly degree variances (see eq. 
(2.18)) scaled to reflect the variance of the thinned and/or merged point 
data with the mean removed was used. Some test computations were 
carried out using empirical covariance functions but this method is 
computationally extremely expensive, due to the large number of data 
involved (see Table 13), and thus was not considered for the final 
predictions. 

b) Bierhammar method (ibid t p. 15) : The 2'x2' Ag were computed as the 

point Ag at the centers of the 2'x2' cells based on the five closest Ag 
to the Ag to be predicted as follows: 

Sj = i if 1 / I ji (5.12) 

where: 1 1 is the distance between the point to be predicted and the 
data point i (power of the prediction = 1). The predictions were also 
carried out using the inverse of the square of the distance as a weight 
(power of prediction = 2): 


*9 = J %/ l — 
3 1=1 *t / 1 = 1 




(5.13) 


LRTITUDE 
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Figure 10. Gravity Predictions in the Vicinity of Hawaii Using the Five 
Closest Points in A Collocation Method. The Grid Spacing is 
2'x2' and the Contour Interval 25 mgal. 
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Figure 11. Gravity Predictions in the Vicinity of Hawaii Using the Five 
Closest Points and the Inverse of the Distance as Weight. 
The Grid Spacing is 2'x2' and the Contour Interval 25 mgal. 
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Figure 12. Gravity Predictions in the Vicinity of Hawaii Using the Five 
Closest Points and the Inverse of the Square of the Distance 
as Weight. The Grid Spacing is 2'x2' and the Contour 
Interval 25 mgal. 
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The following three Figures (10-12) illustrate predictions using the 
three different formulas discussed above: (5.11), (5.12) and (5.13) 
respectively, at the laser stations 7120, 7210 in the vicinity of Hawaii. 
The contour interval is 25 mgal. 


On Table 14 we can see the statistics of the differences between the 
2'x2' mean anomalies predicted using the three different methods in the 
above region. 


Table 14. Statistics of the Differences Between the 2'x2' Mean Anomalies 
Predicted Using: Equation (5.11) (A), Equation (5.12) (B) and 
Equation (5.13) (C) for the Vicinity of Hawaii. _(See Figures 
10, 11, 12). Units are in mgal; Number of 2'x2' Ag Compared: 
22500. 


STATISTICS 

C-B 

C-A 

B-A 

Mean Difference 

- 0.46 

-0.05 

0.41 

RMS Difference 

2.38 

0.97 

2.18 

Standard Deviation Difference 

2.34 

0.97 

2.14 

Maximum Difference 

18.41 

8.32 

10.63 

Minimum Difference 

-10.63 

-4.26 

-18.41 


From the comparisons of Table 14 we can conclude that there is no basic 
difference between the three different prediction methods. Finally, 
method b) with power of prediction = 2 (equation (5.13)) was chosen. 


5.3 Gravity Predictions in the Bermuda Area 

For the transformation of the 15731 available point gravity anomalies 
surrounding the laser station 7067 in Bermuda island (Figure 13) a 
different procedure than the one already described for the data 
surrounding the other laser stations was followed. This procedure used 


is as follows: 











69 


First* the 0*.125 x 0*.125 altimeter anomalies were used to fill-in the 
empty area bounded by the limits 33*20'<4<35* and 30*<4*31*15' in 
latitude and 292*<A<294* and 296*20' <A<298* in longitude. Then, the 
(unthinned) terrestrial point anomalies within the window 
31*15' <f<33*20' and 294*<A<296*20' were used together with the 0*.125 x 
0*.125 anomalies to produce a merged data set. Applying equation (5.13) 
to the merged data set we computed 2'x2' anomalies for the whole 
region (30*<4*35*, 292*<A<298*) except for the central window bounded 
by the limits 32*<4<32*.6 in latitude and 294*.8<A<295*.6 in longitude, 
since the gap that exists within this window (see Figure 13) would give 
unreliable estimates of the 2'x2' anomalies through equation (5.13). 

For this central window the 2'x2' anomalies were computed from the- 
terrestrial point anomalies only if their distance to the five closest 
points ij was less than 3' for every i=l,2,3,4,5. If any of the distances 
*i, i=l,2,3,4,5 was greater than or equal to 3' then the 2'x2' anomaly 
was taken from a 2'x2' anomaly data set computed from 10360 Ohio State 
adjusted GEOS-3/SEASAT altimeter sea surface heights (Liang, 1983) as 
described in Section 8.4. A total of 35 2'x2' anomalies filled the above 
gap. 


The 2'x2' anomaly data set computed as above was then used for 
the computation of the geoid undulation of the laser station 7067. 
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5.4 Terrain Corrections-Indirect Effect 

As it already has been mentioned in Section 2.1, the derivation of 
the Stokes’ equation (equation (2.1)) assumes that no masses external to 
the geoid exist. This is clearly not true for the continents, and thus 
the masses above the geoid have to either be removed from or shifted 
to the geoid. This will result in certain corrections to be applied to the 
computed undulations, known as terrain corrections. For reasons that 
are discussed in (Heiskannen and Moritz, 1967, pp. 151-152) one of the 
most advantageous methods to account for these terrain correction 
computations is the Helmert’s second method of condensation (ibid, p. 
145): The topography is condensed so as to form a surface layer on the 

geoid. The masses are shifted along the local vertical and the total 
mass (of the earth) remains numerically unchanged. This produces an 
attraction change to the free-air anomalies, TC, a potential change V A -V C 
of the topography which results to an indirect effect on the undulation, 
5Nj, and a secondary indirect effect on the gravity anomalies (0.3086 • 
(Wichiencharoen, 1982b). Assuming that the terrain corrections that 
were given for the point anomalies of the continental United States were 
computed using the above Helmert’s second method of condensation, (i.e. 
they equal the TC term discussed above), the terrain-corrections and 
the indirect effects were computed for the continental western United 
States laser stations as follows: To the point free-air anomalies Aj of 

the thinned data sets (see Section 5.1) corrected due to the change of 
the reference ellipsoid (eq. (5.7)) and the effect of the atmosphere (eq. 
(5.9)), the terrain correction that was available was added: 

Sjic = + TC 


(5.14) 
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where Aj TC is the terrain-corrected anomaly 

Aj is the anomaly of the thinned data sets 
TC is the available terrain correction 
Then using the thinned data set (corrected for terrain effects through 
equation (5.14)) the 2'x2' predictions were carried out as described in 
Section 5.2. This method has the disadvantage that for the missing 
2'x2' values no rigorous terrain correction is computed but the terrain 
correction and the value of the free-air anomaly are interpolated as a 
whole from the 5 terrain-corrected anomalies closest to the missing 
value. Then to the computed undulations using the terrain-correctd 
gravity anomalies from equation (5.14) the corresponding indirect effect 
has to be added. The various models for the indirect effect computation 
are discussed in (ibid, p. 19). The simplest and computationally fastest 
but certainly the less accurate model is the model represented by the 
Grushinsky’s formula (ibid, eq. (44)): 

<5N X = = — P — (5.15) 

7 7 

where fiNj: is the indirect effect on the undulation 
V A : potential of the actual topography 
V c : potential of the condensed topography 
G : the Newton’s gravitational constant 

p : the assumed density of the topographic masses (= 2.67 
gr/cm 3 ) 

H : the mean elevation of the l*xl* region surrounding the 

laser station 

7 : the mean value of normal gravity (= 979800 mgal) 
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Table 15 shows the number or the region name of the laser stations 
of the western United States; the maximum terrain correction and the 
arithmetic mean of the terrain correction that was computed from the 
thinned data sets; the mean elevation H of the l*xl* surrounding the 
laser station; and the indirect effect to be added to equation (2.14) if 
the anomalies Aj to be used in the first integral term of (2.14) are taken 
from (5.14). 


Table 15. Information Related to the Terrain Correction and Indirect 
Effect for the Western United States Laser Stations. 


NAME OR 
NUMBER 

TERRAIN CORRECTION (mgal) 


max 

mean 


Wl'lf 

QUI 

81.5 

4.33 

1673 

-0.16 

SAN 

75.4 

2.64 

729 

-0.03 

7082 

57.5 

3.72 

2012 

-0.23 

OVR 

81.5 

5.19 

2363 

-0.32 

GOL 

75.4 

4.97 

674 

-0.03 

FTD 

59.4 

0.99 

1410 

-0.11 

7112 

82.1 

3.91 

1520 

-0.13 

7887 

75.7 

3.49 

100 

-0.00 

HOP 

85.8 

1.63 

1411 

-0.11 

7894 

74.7 

1.97 

122 

-0.00 


The secondary indirect effect on the anomalies (= 0.3086 -fiNi) was 

neglected. 







CHAPTER VI 


ELLIPSOIDAL CORRECTIONS 


Our derivations for the ellipsoidal corrections for all the four 
methods given here follow similar lines as in (Rapp, 1981a). A general 
equation is given for all the four methods, and the ellipsoidal . 
corrections for the traditional methods of Stokes’ and Meissl’s 
modification result as two special cases of the above general equation, 
as it will be shown later in this Chapter. 

The general equation (2.3) that is equivalent to the original Stokes’ 
equation (2.1) is the geoid solution of the boundary value problem and 
is in spherical approximation. Thus the gravity anomaly Ag in (2.3) 
should be substituted with the spherical part of the anomaly, Ag°, where 
(Moritz, 1980, equation (49-21)): 


Ag = Ag° + e 2 Ag x 

(6.1) 

with Ag° = l Ag£ 

(6.2) 

n = o 


Ag 1 = Z Ag x 

(6.3) 

n = o 


Ag® = TT l (A nn) cosmA + B nm sixunX)P nm (sin4>) 
K m=o 

(6.4) 

1 n 

— p X (^nmCOsmA + H nm sinmA ) P nm ( sin+ ) 

m=0 

(6.5) 

ty ■ k -" + x ~ (»;;) + 

(6.6) 
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3(n-3) (n-m-1) (n-m) 
2(2n-3)(2n-l) 

n 3 -3m 2 n-9n 2 -6m 2 -10n+9 
3 (2n+3) (2n-l) 

(3n+5) (n+m+2) (n+m+1) 
2(2n+5) (2n+3) 



_ n ( n-m+1 ) ( n-m+2 ) 
"" = (2n+l)(2n+3) 

n 2 -3m 2 +n 
= (2n+3) (2n—l) 

(n+1 ) ( n+m) ( n+m-1 ) 
nm " (2n+l) (2n-l) 


(6.7) 


( 6 . 8 ) 

(6.9) 


( 6 . 10 ) 


The ’’spherical” undulation from equation (2.3) is (ibid, equation 
49-22): 

N? = || S^cos^Ag’d* + || W 1 (cosy)Ag°d<r (6.11) 

The ellipsoidal form of the equation (2.3) is then (ibid, equation 


(49-26)): 

Nf = N? + AN (6.12) 

with AN = e 2 |^ - ^ sin 2 4]N 0 (6.13) 


N 0 is an approximate value of the undulation which can be taken to be 
the same for all methods, and ♦ is the geodetic latitude of the 
computation point. 


Equation (6.12) can be rewritten using (6.11) and (6.1) as 
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Nf = ^_ }j SiCcosVOUg-e’A gMd<r + || w, (cosV)Ag°d<r + AN 


= 2p- }} s< (cosY0 (Ag-e 2 Ag 1 )d<r + || S { (cos^)Ag°d«r 

^ <r c ny <t-<* c 

+ JJ W,(cosV)Ag°d«r + AN 

7T7 ^ 

= || Sj (cosV)Agd<r + || Sj (cos^)Ag°d<x + W t (cos^) Agd<r 

717 <r c "V-<r c 


+ AN - 1^- || S< (cos^)Ag l d<r (6.14) 

«c 

The first integral term of equation (6.14) is to be computed from the 
terrestrial gravity data and the next three integral terms and 6N will be 
computed from the given set of potential coefficients. 


Thus (6.14) can be rewritten using equation (2.5): 

Nf = 4ny 11 S 1 (cosf)Ajdtf + || S, (cos^)Ag°dff 

<T C O' 

+ 4 ^ || W, (cos^)Ag°d<r + AN - ~ || S 1 (cos^Ag^a 
<r <r c 

= 4^7 11 Si (oosV')AT d <r + ^ J o Q ln Ag° + 2 ; J 0 W >" Ag " + 6ii 

- §“ J1 S^cos^AgM* (6.15) 

Choosing n < M (see Chapter 3), equation (6.15) can be written: 

N ? = 4 ^ 11 S^cos^Ajda + ^ J o (Q in+ W ln )Ag® + ^ J^Q^AgS + AN 
o c 

- J1 S< (cosV)Ag ld <y 

*c 


(6.16) 
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Note that the harmonics Ag£ given by (6.4) will correspond exactly to 
Agg of equation (2.14) if the geocentric latitude ? instead of the geodetic 
latitude ♦ will be used as an argument for the associated Legendre 
functions (Rapp, 1981a, equations (25, (26)). In this case, the 

coefficients 



and (6.16) can be written (ii < M) 

N f = 4 H; ff S,(cos^)ATd<r + ^ J o (Q,„+W in )(n-l)(S) n 
n 

• I (C nm cosmA+S nm sininX ) P nm ( s in$ ) 

m-o 

GM M fal n n 

+ 2 TZ Qi„(n-l)l-J I (C nm cosmA+S nm sinmA)P nm (sin$) 

4,1 7 n = n+l vry m=o 

+ || S,(cosV)Ag *d<r 


(6.17) 


(6.18) 


or 


Nf = Nf + AN + AN, (6.19) 

where Nf is the spherical part of the undulation computed for each of 
the four methods as already discussed in Chapter 2. Also: 

AN, = ~ 11 McosV)Ag»d <r (6.20) 

717 <*c 

We now evaluate the integral of equation (6.20) 

«n, = ' IS 11 Mcos?)A*»*r 

= _ [|| s, (cosViAgMff - || S, (cosV')Ag 1 d<r] 

* <r-<r c 
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= “ [|| S 1 (cosVOAgM* - || S< (cosi>)Ag*d<r] 

717 <r <r 

= - [ 2w ! X «nA«A - 2n ! Q 1n AgA] 

^irry t n =o n=o J 

Dp 2 M 

= - Z (X 1n -Q, n )Ag‘ 

^7 n=o 

= f“ Z (Q 1n -X, n ) J (G nm cosmX+ H nin sinmA)P nm (sin$) (6.21) 

6 7 n = o m-o 

Using equation (6.19), (6.20) and (6.21) we can write equation (6.18) 
as: 

Nf = 4^ JJ s i(cosV)Ag T d«r | o (Q,„+W 1n ) (n-1) (f)" 

n 

• Z (C nm cosmX+S nro sinmA)P nra (sin$) 

m— o 

qm M fal n n 

+ 2^ E. Q in(n-l)|-J Z (C nm cosmA+S nm sirnnX)P nm (sin$) 

“*/ n=n+i m=o 

+ f“ Z (Q 1 n -x in ) z (G nn ,cosmA+H nm sinmX)P nm (sin$) 

n=o m-o 

+ e a (^-f sin a ^N 0 (6.22) 

Equation (6.22) is the desired equation of the geoid undulation 
computation taking into account the ellipsoidal corrections for all the 
four methods. Notice that the fourth correction term (ANj) of equation 
(6.22) arises from the fact that the ellipsoidal effect e a Ag‘ has to be 

A 

removed from the terrestrial anomalies Ag T , whereas the fifth correction 
term (AN) arises from the ellipsoidal correction e a (l/4-3/4sin 2 4)N 0 to the 
spherical undulation N? itself. Taking i=l (Stokes’ method) and starting 
the summation of the second term of the right-hand side in equation 
(6.22) from 2, we take 
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i=1 i . 6. . 22 ^ jj| _ R || (cos^) Ag T d<r 
GM n fa1 n n 

+ 2r> ^ ( Q in + Wi n )(n-l[ 7 J Z (C n(n cosmA+S nm sinmX)P nm (sin$) 

6X 7 n = 2 vr; m— O 

GM * fal n n 

+ 2^Z Q in(n~l)H Z (C„ m cosrnA+S nm sininA)P nm (sin$) 

£tL r n-n+i vrv m=o 

e 2 M n 

+ pZ z (Q xn -X ln ) l ( G„ m cosmA+H n m s inmX ) P n m ( s in$ ) 

+ e 2 ( 4-4 sin 2 *jN 0 (6.23) 

Taking into account the relationships (2.35), (2.36) and (2.39), we have 
Nf = 4 “ }f S(cosV)Ag T d<r 

7 °c 

GM M fain n 

+ 2 7Z Z Q ln (n-l)l-j Z (C nm cosrnA+S nm sinraA)P nm (sin$) 

^ X 7 n— o m=o 

GM H n 

+ nZ z (Qin-Xin) Z (G n(n cosmA+H nm sininA)P n , n (sin$) 

n-o m=o 

+ e * (| - | sin 2 4>)N 0 (6.24) 

Equation (6.24) is exactly the same as equation (31) derived in (ibid, p. 
10845). 

Furthermore for i=2 (Meissl’s method) and starting the summation of 
the second term of the right-hand side in equation (6.22) from 2, we 
take: 

i=2 NS = 4~ |j S 2 (cosf)AgTd<r 

+ 2r? ? ( Q 2 n + W 2 n )(n- 1 ) (f) Z ( C nm cosmA+S nm sinmA )P nm (sin^) 

/ n-2 A m-o 

GM * fal n n 

+ 2r7 n =n+i Q 2 n ^ n_ 1 H7J m | 0 ( c nmCosmA+S nn ,sininA)P nni (sin$) 
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+ f“ I (Q 2 n“ x 2n) l ( Gnm cosmX+ H nn ,sininA )P nm (sin$) 

^7 n=o m=o 

+ e a (| - | sin 2 +)N„ (6.25) 

Taking into account the relationships (2.47), (2.53), (2.57) and (2.59) we 
have: 

HI = \\ (S(cosV)-S 0 )AgTd ff 

flM H f n n 

+ 2^ J |— 1 I (C„ m cosmA+S nin sininX)P nin (sin$) 

^ r 7 n =a m=0 

e 2 M n 

+ o- X (Q 2n -X 2n ) I (G nm cosm\+H nro sinmX)P nni (sin$) 

^7 n~0 m-0 

+ e 2 (| “ | sin 2 *)N 0 (6.26) 

Equation (6.26) is exactly the same as equation (38) derived in (ibid, p. 
10846). 

The equations for the new methods follow in a similar way for i=3 
(Molodenskii’s method) and i=4 (modified Sjoberg’s method) from the 
general equation (6.22). 

To numerically compute the ellipsoidal corrections the OSU86F 
potential coefficient set (C nm , S nm ) up to degree M=36 was used. The 
numerical evaluation of the term 5N was done using an approximation of 
N 0 for all methods taken from 

CM H (aln n 

N ° = ~ I Id l ( C nmCOsmA+S nro s inmX ) P nm ( s in$ ) 

r 7 n= 2 vr/ m=o 


(6.27) 


& 


The numerical evaluation of the term <5Nj was done as follows: First, 
the coefficients {A nm , B nm } were computed using equation (6.8) with {K nm , 
L n(n } = {0, 0} (equation (6.17). Then, the coefficients {G nm , H nm } were 
computed from (6.6) and (6.7). Finally, the term 6N 1 was computed from 
equation (6.21) for all the four methods using the Fourier coefficients 


X < n , Q»„, W ln of each method that have been derived in Chapter 2. 



CHAPTER VII 


LOCAL AVERAGE CORRECTION 


a 

The terrestrial mean anomalies Ag T that are used inside the integral 

A 

of equation (2.14) (Ag T + TC in case that terrain corrections have been 
computed, Ag T otherwise), can be expressed as: 

Ag T = Ag s + fi( Ag 8 ) (7.1) 

where Ag s are the mean anomalies implied by the potential 
coefficient set, used up to degree M 
6( Ag s ) are the errors of the mean anomalies due to the omission 
and commission errors of the given potential coefficients. 

Assuming that the potential coefficients are errorless, the error 
<5(Ag s ) will be due only to the omission errors (truncation of the given 
field at degree M). Taking the average over a spherical cap <r c of both 
sides of eq. (7.1), we obtain: 

| \\ AgTd<r = | JJ Ag s d<r + | JJ<5(Ag 8 )d<r (7.2) 

<*c ff c 

For a high degree field (M = 180, or M = 360) the average of the 
omission errors over a cap of radius on the order of 2* or larger can 
be realistically considered zero. This can be justified by the fact that 
the average of a short wavelength signal (as the omission error 6( Ag s ) 
for a high degree field is) over an area of much longer wavelength 
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(capsize of 2 “) is expected to be zero. 

| Jj fi(Ag*)d<r = 0 (7.3) 

and from (7.2) and (7.3) we must have: 

| || AgT d <r - | |J Ag 8 d<y = 0 (7.4) 

Condition (7.4) expresses the fact that the mean over the cap of 
the terrestrial anomalies should equal to the mean over the cap of the 
anomalies implied by the potential coefficient set. Of course, if the mean 
was taken over the whole earth, then, ideally, both terms of the 
left-hand side of (7.3) would be zero and condition (7.4) would hold 
(Rapp, 1975): 



Ag T d<r - ^ Ag 8 dff = 0-0 = 0 

a 


(7.5) 


In reality 


| |J Ag T d<r - | | J Ag s d<r = AgJ -Agg * 0 
with AgJ = ^ || Ag T d<r 

and Agg = ^ || Ag p d<r 

<*c 


(7.6) 

(7.7) 

(7.8) 


For a numerical verification of (7.6) see columns 7 and 8 of Table 13, 
where 

x * AgJ and x p a Agg 


If we add to the terrestrial anomalies Ag T the difference Agg - AgJ = d, 
then condition (7.4) is met: 
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| Jf (AgT+d)d* - | JjAgsd* 

a c a c 

= S 11 (Ag T +Agg~AgJ)d<r “ ^ 11 Ag s d<r 

*c »c 

= S 11 ^ gTd<r + i 11 dff ~ ^ 11 d«r - ^ |1 Ag B d<r 

tf c a c a c 

= AgJ + Agg - AgJ - Agg = 0 

The difference d = Agg - AgJ is a constant correction to be added to 

a 

the terrestrial anomalies Ag T used to compute the cap contribution. The 
corresponding effect on the geoid undulation can then analytically be 
computed as follows: 


From equation (2.14) with Ag T + d instead of Ag T and n < M we 
obtain: 

N« = 4 “ 11 S 1 (cosy)(AgT+d)d < r + ^ j ( Q jn+ w ln )Ag B 

ff c " ° 


_R 

2 r 


Z_ 

n=o+t 


= 4^ 11 s i( cos V)Ag T dff + g" n | o (Qin + W ln )Agg 

a c " 


*Ty „L, Q '" ie5 + Jl s ' (CO! ” , > dd ’ 

= N, + AN? (7.9) 

where N ( is the geoid undulation computed using each of the four 
methods as discussed in Chapter 2. The local average correction term 


is: 
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dN? = || Sj(cosy)dda 

1X7 <*c 

We have: 

6N? = d || SjCcos^dff 

717 *c 

= d [|| S,(cosy)da - || S< (cosy) da] 
ny a <r-o c 

= d [|| St(cosy)da - || St(cosy)da] (7-10) 

ny <r a 

From equation (2.11), for n=0 we take: 

27iQ j0 = || S f ( cosy) P 0 ( cosy) da, or 
a 

|| S^cosyjda = 27 iQ {0 (7.11) 

a 

Similarly, since X in are the Fourier coefficients of S 1 (cosy) , 

27tXj o = || (cosy) P 0 ( cosy )d<r, or 

a 

|| S i (cosy) da = 27 tX 10 (7.12) 

a 

Substituting (7.11) and (7.12) into (7.10), we obtain: 

<5N? = ^ d (X io -Q io ), or 

6N$ = 2^ (AgS-AgJ)(X io -Q 10 ) (7-13) 

Equation (7.13) is the general expression for the local average 
correction, and it can also be used to analytically compute the effect on 
the undulation due to a constant correction to the gravity data Ag T (e.g. 
a forgotten atmospheric correction). 

For each method, equation (7.13) can be reduced as follows: 
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For i=l (Stokes' method): 

6N? = ~ (Agg-AgJ)(X 10 -Q 10 ) 

= " 2? (AgS-AgJ)Qio (7.14) 

For i=2 (Meissl’s method) 

6N5 = ^ ((Agg-AgX)(X 20 -Q 2 o) 

= ^ (Agg-AgJ) (-2S 0 -Q 2 o) 

= ^ (Agg-AgJ) ( - W 20 -Q 2 o) 

= - 2 ® (Agg-AgJ) (Q 20 +W 20 ) (7.15) 

For i=3 (Molodenskii’s method) 
m = 2 ; (Agg-AgJ) (X 30 -Q 30 ) 

= ^ (Agg-AgJ) (~W 30 ) 

= ~ 2 ; ( Agg-AgJ)W 30 (7.16) 

And for i=4 (Modified Sjoberg’s method) 

5N2 = ^ (Agg-AgJ) (X 40 -Q 4 o) 

= ^ (Agg-AgJ) (-W 40 -Q 40 ) 

= - 2 “ (A«g-AgJ)(W 40 +Q4o) (7.17) 

The numerical values of Q 10 , Q 20 +W 20 » W 30 and W 40 +Q 40 for i>= 2*, n=20, 
M=180, error models (3.1) and (3.2) for the potential coefficients and the 
terrestrial data respectively and model (2.18) for the anomaly degree 
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variances, are: 


Q 10 = -0.075620 
Q2o"^"^2o — 0.035852 

W 30 = -0.053435 
Q 40 +W 40 = -0.0583024 

Then, the corresponding values a t = (R/2r)(Q{ n +W ( n ) in m/mgal of the 

corrections <5N? are (R = 6371000 m and y - 979800 mgal). 

a, = 0.245854 m/mgal 
a 2 = 0.116561 m/mgal 
a 3 = 0.173727 m/mgal 
a 4 = 0.189550 m/mgal 

We can see that Meissl’s method is the least sensitive and Stokes’ 

method is the most sensitive to a constant correction. Molodenskii’s and 

modified Sjoberg’s methods have sensitivities that are between the 

Stokes’ and Meissl’s methods. For the numerical evaluation of <5N?, the 

values of Ag§, AgJ were taken from Table 13. 


(7.19) 


(7.18) 


C -9- 



CHAPTER VIII 


NUMERICAL RESULTS AND STATISTICS 
8.1 Computational details 

As discussed in Chapter 5 the numerical integration of the functions 
S f (cos i>) was computed using 2'x2', 6'xlO' and 10'xlO' mean terrestrial 
Ag T . The selected capsize was 2* and the functions Sj(cos i>) were first 
tabulated at a spherical interval of 1" out to a spherical distance of 3*. 
Then the evaluation of S t (cosf) for each block was done by numerical 
integration as follows: Each block was subdivided into elements 

depending on the separation f between the computation point (laser 
station) and the cell. For tp- 4' the number of subdivisions was 64; for 
4'<^<8' we had 16 subdivisions; for 8'<^12' the subdivisions were 4 
and for 12'<Y>^2* the evaluation was made at the center point of the cell 
(Engelis et al., 1985a). The general equation including ellipsoidal 
corrections (eq. (6.22)) was applied for all the methods. The outer zone 
contribution (second and third term in (6.22)) was computed from the 
OSU86F potential coefficients taken up to degree M=180. Since there is 
no zero or first degree potential coefficients in the OSU86F field the 
summations started from 2, although the theoretical formulas developed 
start their summations from zero. The correct procedure would be to 
start our summations from zero. The ellipsoidal corrections were 
computed up to degree M=36 as mentioned in Chapter 6. The models for 
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the error degree variances of the terrestrial anomalies and the anomalies 
implied by the potential coefficients are important for the error analysis 
of the four methods, but they also play a very important role for the 
modified Sjoberg’s method, since for this method the Fourier coefficients 
X 4n , W 4n and Q 4n are based on them. The model for the anomaly error 
degree variances of the potential coefficient set was taken from equation 
(3.1). The models for the terrestrial anomaly error degree variances 
were selected for each region as Table 16 shows: 

Table 16. Models for the Terrestrial Error Anomaly Degree Variances 
Selected by Region. 


Region 

Model 

UNITED STATES 

BERMUDA 

AUSTRALIA 

EUROPE 

JAPAN 

B.l 

B.l 

B.l 

A [0 n =0 n (6'xlO')] 

A r0„=0„aO'xlO')l 


Finally the number n of harmonics that were removed from the Stokes’ 
function for the Molodenskii’s and the modified Sjoberg’s method was 
taken to be 20 (<M=180). For comparison purposes the point undulation 
of the laser station was computed using the OSU86F potential 
coefficients in equation (6.23) with M=360. 

8.2 Numerical results 

Table 17 shows the results of the undulation computations with the 
Stokes’ method. The first and second columns show the number and the 
name of the laser stations; the third column shows the reference for the 
ellipsoidal height of the stations; the next three columns represent the 
three components of the undulation (see equation (6.22)); cap 

contribution, NC, potential coefficient contribution (outer zone) NP, and 
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Table 17. Geoid Undulations of the Laser Stations Using Stokes’ Method 
(Units are in meters). 


NLM 

NAME 

REF 

NC 

NP 

NE 

N 

OSU86 F 

ELL-ORT 

7051 

0U I LAS 

BM 

2.79 

-27.12 

-0.01 

-24.35 

-23.64 

-22.57 

7109 

QUILAS 

BM 

2.79 

-27,13 

-0.01 

-24.35 

-23.65 

-22.59 

7836 

OUILAS 

UN 

2.79 

-27.13 

-0.01 

-24.35 

-23.65 

-19.30 

7062 

SANOIE 

BM 

-1.41 

-33. 88 

-0.01 

-35.30 

-34.22 

-32.47 

71 10 

MONLAS 

8M 

-0.33 

-33.06 

-0.01 

-33.39 

-32.78 

-30.86 

7082 

BEARLK 

BM 

1.75 

-18.74 

0.00 

-16.99 

-15.16 

-12.58 

7084 

OVRLAS 

BM 

1.22 

-28.63 

-0.01 

-27 . 42 

-26-29 

-24.90 

7114 

OVRLAS 

8M 

1.21 

-28.62 

-0.01 

-27.42 

-26.30 

-24.84 

7085 

GOLOLS 

BM 

-1.76 

-30. 17 

-0.01 

-31.94 

-31.70 

-29.68 

7115 

GOLLAS 

BM 

-1.85 

-30.31 

-0.01 

-32.16 

-31.84 

-29.58 

7036 

FTOAVS 

BM 

1.51 

-23.69 

0.01 

-22.17 

-21.36 

-20.78 

7835 

FORLAS 

BM 

1.51 

-23.69 

0.01 

-22.17 

-21.36 

-20.83 

7112 

PLALAS 

BM 

2.40 

-21.37 

0.00 

-18 . 96 

-18.28 

-17.40 

7837 

VANLAS 

UN 

-2.68 

-33.43 

-0,01 

-36 . 12 

-36.47 

-31.14 

7838 

HOPLAS 

UN 

-2.01 

-28 . 46 

0.00 

-30.48 

-29.48 

-2.20 

7921 

HOPLAS 

UN 

-2.02 

-28 . 46 

0.00 

-30. 48 

-29.48 

-29.47 

7894 

YUMLAS 

UN 

-4.34 

-30.00 

-0.01 

-34.34 

-33.56 

-28.63 

7063 

STALAS 

BM 

0.04 

-33.26 

-0.01 

-33.23 

-33.08 

-32.92 

7064 

GSFCLS 

BM 

0.04 

-33.26 

-0,01 

-33.23 

-33.08 

-33.64 

7 100 

GSF100 

BM 

0.02 

-33.26 

-0.01 

-33.25 

-33.09 

-32.50 

7101 

GSF 101 

BM 

0.02 

-33.26 

-0.01 

-33.25 

-33.09 

-33.08 

7102 

GSF 102 

BM 

0.04 

-33.26 

-0.01 

-33.23 

-33.08 

-32 . 89 

7 103 

GSF 103 

BM 

0.04 

-33.26 

-0.01 

-33.23 

-33.08 

-32.93 

7104 

GSF 104 

BM 

0.03 

-33.26 

-0.01 

-33.25 

-33.08 

-32.94 

7105 

GSF 105 

BM 

0.04 

-33.26 

-0.01 

-33.23 

-33.08 

-32.93 

7069 

RAMLAS 

BM 

1.02 

-30.25 

0.01 

-29.22 

-29.52 

-29 . 13 

709 1 

HAYLAS 

BM 

0.54 

-29.54 

-0.02 

-29.02 

-28.61 

-27.55 

7120 

MU I LAS 

BM 

20.7 7 

-1.83 

0.00 

18.93 

13.02 

20.46 

7210 

MAULAS 

IA 

20.80 

-1.83 

0.00 

18.96 

13.02 

20.4 9 

7090 

YARLAS 

BM 

-1.00 

-23.77 

0.00 

-24.77 

-25.02 

-24.28 

7943 

ORRLAS 

UN 

7.47 

1 1 .08 

0.02 

18.57 

20.36 

20.38 

7067 

B 0 I LAS 

BM 

5.69 

-37 . 73 

-0.01 

-32.05 

-35.23 

-31.81 

7805 

FINLAS 

IA 

-2.54 

23.07 

-0.01 

20.53 

20.54 

19.95 

7939 

MATLAS 

UN 

9.91 

36.63 

0.01 

46.55 

45.66 

46.39 

7035 

GRALAS 

UN 

6 . 98 

44.93 

0.02 

51.93 

51 .04 

42 . 85 

7840 

RGOLAS 

UN 

-2.46 

48.19 

0.01 

45.75 

4 6.06 

45.58 

0833 

K00LAS 

UN 

-0.02 

44.14 

0.01 

44 . 13 

44 . 80 

44 .66 

7834 

WETLAS 

IA 

4.49 

42.17 

0.02 

46.68 

47.02 

47.65 

7338 

SHOLAS 

UN 

6.85 

32.05 

0.03 

38.93 

37.60 

39.97 


Table 18. Geoid Undulations of the Laser Stations Using Meissl’s Method 
(Units are in meters). 


NLM 

NAME 

REF 

NC 

NP 

NE 

N 

OSU06 F 

ELL-ORT 

7051 

QUILAS 

BM 

1.73 

-26 .07 

0.00 

-24.34 

-23.64 

-22.57 

7109 

QUILAS 

BM 

1 .73 

-26.07 

0.00 

-24.34 

-23.65 

-22.59 

7886 

QUILAS 

UN 

1.73 

-26.07 

0.00 

-24.34 

-23.65 

-19.30 

7062 

SANOIE 

BM 

-0.05 

-34.85 

0.00 

-34.90 

-34.22 

-32 . 47 

7110 

MONLAS 

BM 

1.13 

-34.05 

0.00 

-32.87 

-32.78 

-30.86 

7032 

BEARLK 

BM 

0.08 

-16.44 

0.00 

-16.36 

-15.16 

-12.58 

7084 

OVRLAS 

BM 

0.92 

-28.50 

0.00 

-27.59 

-26.29 

-24.90 

7114 

OVRLAS 

BM 

0.91 

-28.50 

0,00 

-27.59 

-26.30 

-24 . 84 

7085 

GOLDLS 

BM 

-1.31 

-30.50 

0.00 

-31 .81 

-31.70 

-29.68 

7115 

GOLLAS 

BM 

-1.25 

-30.72 

0.00 

-31.97 

-31.84 

-29.58 

7086 

FTOAVS 

BM 

1 . 59 

-23.23 

0.00 

-21 . 64 

-21 .36 

-20.78 

7885 

FORLAS 

BM 

1.59 

-23.23 

0.00 

-21.64 

-21.36 

-20.83 

7112 

PLALAS 

BM 

-0.18 

-18.60 

0.00 

-18.77 

-18.28 

-17.40 

7887 

VANLAS 

UN 

-0. 73 

-35.46 

0.00 

-36. 19 

-36.47 

-31.14 

7888 

HOPLAS 

UN 

-0.80 

-29 . 26 

0.00 

-30.07 

-29.48 

-2.20 

7921 

HOPLAS 

UN 

-0.81 

-29.26 

0.00 

-30.07 

-29.48 

-29.47 

7894 

YUMLAS 

UN 

-2.14 

-31 .71 

0.00 

-33.86 

-33.56 

-28.63 

7063 

STALAS 

BM 

0.36 

-33.88 

-0.01 

-33.03 

-33.08 

-32 . 92 

7064 

GSFCLS 

8M 

0.36 

-33.88 

-0.01 

-33.03 

-33.08 

-33.64 

7 100 

GSF 100 

BM 

0.84 

-33.88 

-0.01 

-33.05 

-33.09 

-32 . 50 

7101 

GSF 101 

BM 

0.84 

-33.88 

-0.01 

-33.05 

-33.09 

-33.08 

7102 

GSF 102 

BM 

0.86 

-33.88 

-0.01 

-33.03 

-33.08 

-32.89 

7103 

GSF I 03 

BM 

0.86 

-33 .88 

-0.01 

-33.03 

-33.08 

-32.93 

7104 

GSF 104 

BM 

0.85 

-33.88 

-0.01 

-33.04 

-33.08 

-32 . 94 

7105 

GSF 105 

BM 

0.86 

-33.88 

-0.01 

-33.03 

-33.08 

-32 . 93 

7069 

RAMLAS 

BM 

0.50 

-30.04 

0.00 

-29.53 

-29.52 

-29.13 

7091 

HAYLAS 

BM 

0.49 

-29. 17 

-0.01 

-28.71 

-28.61 

-27.55 . 

7120 

MU I LAS 

BM 

14.70 

3.21 

0.00 

17.91 

13.02 

20.46 

7210 

MAULAS 

IA 

14.73 

3.20 

0.00 

17.94 

13.02 

20.49 

7090 

YARLAS 

BM 

-0.8 7 

-24.21 

0.00 

-25.08 

-25.02 

-24 . 28 

7943 

ORRLAS 

UN 

4.44 

14.72 

0.01 

19. 17 

20.36 

20.38 

7067 

BO I LAS 

BM 

6.63 

-38.72 

-0.01 

-32. 10 

-35.23 

-31.81 

7805 

FINLAS 

IA 

-0.53 

21.20 

0.00 

20.67 

20.54 

19.95 

7939 

MATLAS 

UN 

6.25 

40. 12 

0.00 

46.37 

45.66 

46.39 

7835 

GRALAS 

UN 

5.02 

46.81 

0.01 

51.83 

51 .04 

42.85 

7840 

RGOLAS 

UN 

-1.53 

47 . 40 

0.01 

45.88 

46.06 

45.58 

8833 

KOOLAS 

UN 

-0.06 

44.55 

0.01 

44 . 49 

44.80 

44 . 66 

7834 

WETLAS 

IA 

2.28 

44.65 

0.01 

46.94 

47.02 

47.65 

7838 

SHOLAS 

UN 

4.79 

34.38 

0.01 

39. 18 

37.60 

39.97 
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ellispoidal correction NE (tfN 1 +5N terms of eq. (6.22)). For the 17 laser 
stations in the Western United States the NC includes the cap 
contribution using terrain-corrected gravity anomalies and the indirect 
effect 5Nj (see eq. (5.15)): 

NC = T~ II g T +TC)S(cosf)d<r + 5Ni (8.1) 

7X7 a o 

In Tables 17 through 21 the "NC" column will always be given by 

(8.1) for the various modifications of the Stokes’ function S^cos 1 >), 
unless otherwise stated. The seventh colunm (N) is the total undulation 
computed using the OSU86F set up to M-360 (equation (6.27)). Finally, 
the last column is the value of the undulation computed from equation 

(4.2) (Ellipsoidal minus orthometric height). Units are in meters. Tables 
18, 19 and 20 show the results for Meissl’s, Molodenskii’s and modified 
Sjoberg’s methods and their description is identical with the description 
of Table 17 above. 

Table 21 summarizes the results of the four methods by region: The 

first two columns show the number and name of the laser station; 
columns three to six show the total undulation computed using the four 
different methods: Stokes’, Meissl’s, Molodenskii’s and modified Sjoberg’s 
method respectively; the seventh column shows the height reference; the 
eighth column shows the value of the undulation using the OSU86F set 
as above and the last column shows the value of the undulation 
computed from equation (4.2). 


OF POOR 


* i < -i;+.i I • J 

quality 
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Table 19. 


Geoid Undulations of the Laser Stations Using Molodenskii’s 
Method (Units are in meters). 


MJM 

NAME 

REF 

NC 

NP 

NE 

N 

0SU86 F 

ELL-ORT 

7051 

QUILAS 

BM 

2.2* 

-26.53 

-0.01 

-24.34 

-23.64 

-22.57 

7109 

QUILAS 

BM 

2.2* 

-26.53 

-0.01 

-24.34 

-23.65 

-22.59 

7886 

QUILAS 

UN 

2.2* 

-26.53 

-0.01 

-24.34 

-23.65 

-19.30 

7062 

SAND IE 

BM 

-8.65 

-34.41 

0.00 

-35.07 

-34.22 

-32.47 

7110 

MONLAS 

BM 

0.51 

-33.60 

0.00 

-33. 10 

-32.78 

-30.86 

7032 

BEAR L K 

BM 

0.82 

-17.46 

0.00 

-16.64 

-15.16 

-12.58 

7084 

OVRLAS 

BM 

1 .04 

-28.55 

-0.01 

-27.52 

-26.29 

-24.90 

71 14 

OVRLAS 

BM 

t .03 

-28.55 

-0.01 

-27.52 

-26.30 

-24 . 84 

7895 

GOLDLS 

BM 

-1.51 

-30.36 

-0.01 

-31.87 

-31 . 70 

-29.68 

71 15 

GOLLAS 

BM 

-1.51 

-30.54 

0.00 

-32.06 

-31 .84 

-29.58 

7036 

FTDAVS 

BM 

1.55 

-23.43 

0.01 

-21 . 87 

-21 .36 

-20.78 

7895 

FORLAS 

BM 

1.55 

-23.43 

0.01 

-21.87 

-21 .36 

-20.83 

7112 

PLALAS 

BM 

0.96 

-19.82 

0.00 

-18.86 

-18.28 

-17.40 

7887 

VANLAS 

UN 

-1.59 

-34.56 

-0.01 

-36 . 16 

-36.47 

-31.14 

7838 

HOPLAS 

UN 

-1.34 

-28.91 

0.00 

-30.25 

-29 . 48 

-2.20 

7921 

HOPLAS 

UN 

-1.34 

-28.91 

0.00 

-30.25 

-29 . 48 

-29 .47 

7394 

YUMLAS 

UN 

-3.11 

-30.96 

0.00 

-34.07 

-33.56 

-28'.63 

7063 

STALAS 

BM 

0. 49 

-33.60 

-0.01 

-33.12 

-33.08 

-32.92 

7064 

GSFCLS 

BM 

0. 49 

-33.60 

-0.01 

-33.12 

-33.08 

-33.64 

7100 

GSF 100 

BM 

0.48 

-33.60 

-0.01 

-33. 14 

-33.09 

-32.50 

7101 

GSF101 

BM 

0.48 

-33.60 

-0.01 

-33. 14 

-33.09 

-33.08 

7102 

GSF 102 

8M 

0.49 

-33.60 

-0.01 

-33 . 12 

-33.08 

-32.89 

7103 

GSF 1 03 

BM 

0.49 

-33.60 

-0.01 

-33. 12 

-33.08 

-32 . 93 

7104 

GSF 1 04 

BM 

0.48 

-33.60 

-0.01 

-33 .13 

-33.08 

-32.94 

7105 

GSF 1 05 

BM 

0.49 

-33.60 

-0.01 

-33.12 

-33.08 

-32 . 93 

7069 

RAMLAS 

BM 

0.73 

-30.13 

0.01 

-29 . 39 

-29.52 

-29 . 13 

7091 

HAYLAS 

BM 

0.51 

-29.34 

-0.01 

-28 . 85 

-28.61 

-27.55 

7 170 

MU ILAS 

BM 

17.36 

1 .00 

0.00 

18.36 

13.02 

20.46 

7210 

MAULAS 

IA 

17.40 

1 .00 

0.00 

18.39 

13.02 

20.4 9 

7090 

YARLAS 

BM 

-0.83 

-24.01 

0.00 

-24.84 

-25.02 

-24.23 

7943 

ORRLAS 

UN 

5.77 

13.12 

0.01 

13.91 

20.36 

20.33 

7067 

BO I LAS 

BM 

6.21 

-38.28 

-0.01 

-32.08 

-35.23 

-31.81 

7905 

F INLAS 

1A 

-1.42 

22.03 

-0.01 

20.60 

20.54 

19.95 

7939 

MATLAS 

UN 

7.86 

38.59 

0.00 

46 . 45 

45.66 

46 .39 

78 3 5 

GRALAS 

UN 

5.87 

45.99 

0.01 

51.87 

51 .04 

42 . 85 

7340 

RGOLAS 

UN 

-1.94 

47.75 

0.01 

45.81 

46.06 

45 . 59 

3833 

KOOLAS 

UN 

-0.05 

44.37 

0.00 

44.32 

44 . 80 

44.66 

7834 

WETLAS 

IA 

3.26 

43.55 

0.01 

46 . 82 

47 . 32 

47 . 65 

7838 

SHOLAS 

UN 

5.70 

33.35 

0.02 

39 .07 

37.60 

39 . 97 

Table 20, Geoid 

Undulations 

of 

the Laser 

Stations 

Using Modified 


Sjoberg's Method (Units are in meters). 


nuM 

NAME 

REF 

NC 

NP 

NE 

N 

OSU86 F 

ELL-ORT 

705 1 

QUILAS 

BM 

2.35 

-26 . 60 

-0.01 

-24.26 

-23.64 

-22.57 

7109 

QUILAS 

BM 

2.35 

-26.61 

-0.01 

-24.27 

-23.65 

-22.59 

7836 

QUILAS 

UN 

2.35 

-26.61 

-0.01 

-24.26 

-23.65 

-19.30 

7062 

SANOIE 

BM 

-0.33 

-34.20 

0.00 

-35.04 

-34.22 

-32.47 

71 10 

MONLAS 

BM 

0.32 

-33.37 

0.00 

-33.06 

-32.78 

-30.86 

7082 

B EAR t K 

BM 

1.17 

-17.77 

0.00 

-16.60 

-15 . 16 

-12.58 

7084 

OVRLAS 

BM 

0.36 

-28.43 

-0.01 

-27.57 

-26.29 

-24.90 

7114 

OVRLAS 

BM 

0.85 

-28.42 

-0.01 

-27.58 

-26.30 

-24 . 84 

7035 

GOLDLS 

BM 

- 1.58 

-30.32 

-0.01 

-31.91 

-31.70 

-29.68 

7115 

GOLLAS 

BM 

-1.59 

-30. 49 

0.00 

-32.08 

-31 .84 

-29.58 

7036 

FTDAVS 

BM 

1.54 

-23.45 

0.01 

-21.90 

-21 . 36 

-20.78 

7385 

FORLAS 

BM 

1.54 

-23.45 

0.01 

-21.90 

-21.36 

-20.83 

7112 

PLALAS 

BM 

1.23 

-20. 18 

0.00 

-18.95 

-18.28 

-17.40 

7837 

VANLAS 

UN 

-1.87 

-34 . 26 

-0.01 

-36 . 14 

-36.47 

-31.14 

7338 

HOPLAS 

UN 

-1.54 

-28.78 

0.00 

-30.32 

-29.48 

-2 . 20 

7921 

HOPLAS 

UN 

-1.54 

-28 . 78 

0.00 

-30.32 

-29.48 

-29 .47 

7894 

YUMLAS 

UN 

-3.23 

-30.86 

0.00 

-34 .09 

-33.56 

-28 . 63 

7063 

STALAS 

BM 

0.40 

-33.50 

-0.01 

-33.11 

-33.08 

-32 . 32 

7064 

GSFCLS 

BM 

0.40 

-33 . 50 

-0.01 

-33 . 10 

-33.08 

-33 . 64 

7 100 

GSF 100 

BM 

0.39 

-33.50 

-0.01 

-33.13 

-33.09 

-32.50 

7101 

GSF 101 

BM 

0.39 

-33.50 

-0.01 

-33.12 

-33.09 

-33.08 

7102 

GSF 102 

BM 

0.40 

-33.50 

-0.01 

-33. 10 

-33.08 

-32 . 89 

7103 

GSF 103 

BM 

0.40 

-33.50 

-0.01 

-33. 10 

-33.08 

-32.93 

7104 

GSF 104 

BM 

0.39 

-33.50 

-0.01 

-33.12 

-33.08 

-32.94 

7105 

GSF 1 05 

BM 

0.40 

-33.50 

-0.01 

-33.11 

-33.08 

-32.93 

7063 

RAMLAS 

BM 

0.35 

-30. 18 

0.01 

-29 . 33 

-29.52 

-29. 13 

709 1 

HAYLAS 

BM 

0.45 

-29.36 

-0.01 

-28.92 

-28.61 

-27.55 

7120 

MU I LAS 

BM 

17.77 

0.86 

0.00 

18.63 

13.02 

20.46 

7210 

MAULAS 

IA 

17.30 

0.86 

0.00 

18.66 

13.02 

20.49 

7090 

YARLAS 

BM 

-0.97 

-23.91 

0.00 

-24 . 88 

-25.02 

-24 . 28 

7943 

ORRLAS 

UN 

5.90 

12.94 

0.01 

18.85 

20.36 

20.38 

7067 

BO I LAS 

BM 

5.99 

-37.95 

-0.01 

-31.97 

-35.23 

-31.81 

7805 

F INLAS 

IA 

-1.80 

22.72 

-0.01 

20.91 

20.54 

19 . 95 

7939 

MATLAS 

UN 

7.9 1 

37.84 

0.01 

45 . 76 

45.66 

46.39 

7835 

GRALAS 

UN 

5.90 

45.64 

0.01 

51.55 

51.04 

42.85 

7840 

RGOLAS 

UN 

-1.94 

47.89 

0.01 

45 . 96 

46.06 

45.58 

8333 

KOOLAS 

UN 

0.00 

44.21 

0.01 

44.22 

44.80 

44.66 

7834 

WETLAS 

IA 

3.62 

42.73 

0.02 

46.37 

47.02 

47.65 

7838 

SHOLAS 

UN 

5.99 

32.61 

0.02 

38.63 

37.60 

39.97 
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Table 21. Total Geoid Undulations of the Laser Stations Using All the 
Pour Methods (Units are in meters). 


NUM 

NAME 

STOKES 

ME I SSL 

MOLOO. 

SJOB . 

REF 

OSU86F 

ELL-ORT 




WEST UNITED STATES 




7051 

QUILAS 

-24.35 

-24.34 

-24.34 

-24.26 

BM 

-23.64 

-22.57 

7109 

QUILAS 

-24,35 

-24.34 

-24.34 

-24.27 

BM 

-23,65 

-22.59 

7886 

QUILAS 

-24.35 

-24.34 

-24.34 

-24.26 

UN 

-23.65 

-19.30 

7062 

SAND I E 

-35.30 

-34.90 

-35.07 

-35.04 

BM 

-34.22 

-32.47 

7110 

MONLAS 

-33.39 

-32.87 

-33.10 

-33.06 

BM 

-32.78 

-30.86 

7082 

BEARLK 

-16.99 

-16.36 

-16.64 

-16.60 

BM 

-15. 16 

-12.58 

7084 

OVRLAS 

-27.42 

-27.59 

-27.52 

-27.57 

BM 

-26.29 

-24.90 

7114 

OVRLAS 

-27.42 

-27.59 

-27.52 

-27.58 

BM 

-26.30 

-24.84 

7085 

GOIDLS 

-31.94 

-31.81 

-31.87 

-31.91 

BM 

-31 .70 

-29.68 

7115 

GOLLAS 

-32.16 

-31.97 

-32.06 

-32.08 

BM 

-31 .84 

-29.58 

7086 

FTDAVS 

-22.17 

-21.64 

-21.87 

-21.90 

BM 

-21.36 

-20.78 

7885 

FORLAS 

-22.17 

-21 .64 

-21 .87 

-21.90 

BM 

-21 .36 

-20.83 

7112 

PLALAS 

-18.96 

-19.77 

-18.86 

-18.95 

BM 

-13.20 

-17.40 

7887 

VANLAS 

-36.12 

-36.19 

-36. 16 

-36. 14 

UN 

-36.47 

-31.14 

7888 

HOPLAS 

-30.48 

-30.07 

-30.25 

-30.32 

UN 

-29.48 

-2 . 20 

792 1 

HOPLAS 

-30.48 

-30.07 

-30.25 

-30.32 

UN 

-29.48 

-29.47 

7894 

YUMLAS 

-34.34 

-33.86 

-34.07 

-34 .09 

UN 

-33.56 

-28 . 63 




EAST UNITED STATES 




7063 

STALAS 

-33.23 

-33.03 

-33. 12 

-33. 1 1 

BM 

-33.08 

-32.92 

7064 

GSFCLS 

-33.23 

-33.03 

-33. 12 

-33. 10 

BM 

-33.08 

-33.64 

7100 

GSF100 

-33.25 

-33.05 

-33. 14 

-33. 13 

BM 

-33.09 

-32.50 

7101 

GSF101 

-33.25 

-33.05 

-33. 14 

-33.12 

BM 

-33.09 

-33.08 

7102 

GSF102 

-33.23 

-33.03 

-33.12 

-33.10 

BM 

-33.08 

-32.89 

7103 

GSF103 

-33.23 

-33.03 

-33. 12 

-33. 10 

BM 

-33.08 

-32.93 

7104 

GSF104 

-33.25 , 

-33.04 

-33.13 

-33. 12 

BM 

-33.08 

-32.94 

7105 

GSF105 

-33.23 

-33.03 

-33.12 

-33. I 1 

BM 

-33.08 

-32 . 93 

7069 

RAMLAS 

-29.22 

-29.53 

-29.39 

-29.33 

BM 

-29.52 

-29. 13 

7091 

HAYLAS 

-29.02 

-28.71 

-28.85 

-28.92 

BM 

-28 .61 

-27.55 




HAWA I I 






7120 

MU I LAS 

18.93 

17.91 

18.36 

18.63 

BM 

13.02 

20.46 

7210 

MAULAS 

18.96 

17.94 

18.39 

18.66 

IA 

13.02 

20.49 




AUSTRALIA 





7090 

YARLAS 

-24.77 

-25.08 

-24.84 

-24.88 

BM 

-25.02 

-24.28 

7943 

ORRLAS 

18.57 

19.17 

18.91 

18.85 

UN 

20.36 

20.38 




BERMUDA 






7067 

BD I LAS 

-32.05 

-32.10 

-32.08 

-31.97 

BM 

-35.23 

-31.81 




EUROPE 






7805 

FINLAS 

20.53 

20.67 

20.60 

20.91 

IA 

20.54 

19.95 

7939 

MATLAS 

46.55 

46.37 

46.45 

45.76 

UN 

45.66 

46.39 

7835 

GRALAS 

51.93 

51.83 

51.87 

51.55 

UN 

51 .04 

42.85 

7840 

RGOLAS 

45.75 

45.88 

45.81 

45.96 

UN 

46.06 

45.58 

8833 

KOOLAS 

44.13 

44.49 

44.32 

44.22 

UN 

44 . 80 

44.66 

7834 

WETLAS 

46.68 

46.94 

46.82 

46.37 

IA 

47.02 

47 . 65 




JAPAN 






7038 

SHOLAS 

38.93 

39.18 

39.07 

38.63 

UN 

37.60 

39 .97 
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In Tables 22 and 23 the undulations for the western United States 
laser stations computed using terrain-corrected (NC from (8.1)) and 
uncorrected gravity data with the modified Sjoberg’s method are shown. 
In Table 24 the total geoid undulation of the laser stations using the 
local average correction discussed in Chapter 7 for all the four methods 
are shown. The magnitude of this effect is 1-2 m. It would be helpful 
to have a table of residuals between the undulations computed using all 
the four methods and the undulations computed from equation (4.2). 
These residuals were computed only for the 28 laser stations with known 
height references. Table 25 shows the number and name of the laser 
stations, the differences Nt-fh-H) where are given in Table 21 for all 

the four methods and h-H is the undulation computed from eq. (4.2), the 
height references and the residuals between the undulations computed 
from the OSU86F set (M=360) and the h-H value. 

The corresponding table of residuals when the local average 
correction is included is shown in Table 26. The statistics of the 
residuals between the terrain-corrected Njc and the uncorrected N un 
using the modified Sjoberg’s method and the N T = h-H value are shown 
in Table 27 for the 12 laser stations in the western United States 
(corresponded values in Tables 22 and 23). 

In Table 28 the statistics of the residuals of the Table 25 are shown 
(undulations computed using all four methods vs. the N T value) and in 
Table 29 the corresponded statistics including the local average 


correction are shown. 
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Table 22. Geoid Undulations of the 17 Laser Stations in the Western 
United States Using Modified Sjoberg’s Method and 
Terrain-corrected Gravity Anomalies (Units are in meters). 


NUM 

NAME 

REF 

NC 

NP 

NE 

N 

0SU86 F 

ELL-ORT 

705 1 

QUILAS 

BM 

2.35 

-26.60 

-0.01 

-24.26 

-23.64 

-22.57 

7103 

QUILAS 

BM 

2.35 

-26 .61 

-0.01 

-24.27 

-23.65 

-22.59 

788 S 

QUILAS 

UN 

2.35 

-26.61 

-0.01 

-24.26 

-23.65 

-19.30 

7062 

SAND I E 

BM 

-0.83 

-34 . 20 

0.00 

-35.04 

-34.22 

-32.47 

7110 

MONLAS 

BM 

0.32 

-33.37 

0.00 

-33.06 

-32.78 

-30.86 

7082 

BEARLK 

BM 

1.17 

-17.77 

0.00 

-16.60 

-15. 16 

-12.58 

708 4 

OVRLAS 

BM 

0.86 

-28.43 

-0.01 

-27.57 

-26.29 

-24.90 

7114 

OVRLAS 

BM 

0.85 

-28.42 

-0.01 

-27.58 

-26.30 

-24.84 

7085 

GOLDLS 

BM 

-1.58 

-30.32 

-0.01 

-31.91 

-31.70 

-29.68 

7115 

GOLLAS 

BM 

-1.59 

-30. 49 

0.00 

-32.08 

-31.84 

-29.58 

7086 

FTDAVS 

BM 

1 .54 

-23.45 

0.01 

-21.90 

-21.36 

-20.78 

7885 

FORLAS 

BM 

1.54 

-23.45 

0.01 

-21 .90 

-21.36 

-20.83 

7112 

PLALAS 

BM 

1 .23 

-20. 18 

0.00 

-18.95 

-18.23 

-17.40 

7887 

VANLAS 

UN 

-1 .87 

-34.26 

-0.01 

-36.14 

-36.47 

-31.14 

7883 

HOPLAS 

UN 

-1 .54 

-28.73 

0.00 

-30.32 

-29.48 

-2.20 

7921 

HOPLAS 

UN 

-1.54 

-28.73 

0.00 

-30.32 

-29.48 

-29.47 

7894 

YUMLAS 

UN 

-3.23 

-30.36 

0.00 

-34.09 

-33.56 

-23.63 


Table 23. Geoid Undulations of the 17 Laser Stations in the Western 
United States Using Modified Sjoberg’s Method and 
Uncorrected Gravity Anomalies (Units are in meters). 


NUM 

NAME 

REF 

NC 

NP 

NE 

N 

0SU86 F 

ELL-ORT 

705 1 

QUILAS 

BM 

1 .77 

-26.60 

-0.01 

-24 . 84 

-23.64 

-22.57 

7109 

QUILAS 

BM 

1.77 

-26.61 

-0.01 

-24.84 

-23 . 65 

-22 . 59 

7886 

QUILAS 

UN 

1 .77 

-26.61 

-0.01 

-24 .84 

-23 . 65 

-19.30 

7062 

SAND IE 

BM 

-1.22 

-34.20 

0.00 

-35.43 

-34 . 22 

-32 . 47 

71 10 

MONLAS 

BM 

-0.07 

-33.37 

0.00 

-33.45 

-32 . 78 

-30.86 

7082 

BEARLK 

BM 

0.87 

-17.77 

0.00 

-16.90 

-15. 16 

-12.58 

7084 

OVRLAS 

BM 

-0. 16 

-28.43 

-0.01 

-28.59 

-26 . 29 

-24.90 

7114 

OVRLAS 

BM 

-0.17 

-23.42 

-0.01 

-28.60 

-26 . 30 

-24.84 

7035 

GOLDLS 

BM 

-2.21 

-30.32 

-0.01 

-32.54 

-31.70 

-29 . 68 

7115 

GOLLAS 

BM 

-2.16 

-30.49 

0.00 

-32.65 

-31.84 

-29 . 58 

7086 

FTDAVS 

BM 

1.32 

-23.45 

0.01 

-22.12 

-21.36 

-20.78 

7885 

FORLAS 

BM 

1.32 

-23.45 

0.01 

-22.12 

-21 .36 

-20.83 

7112 

PLALAS 

BM 

1 .05 

-20. 18 

0.00 

-19.13 

-18.28 

-17.40 

7887 

VANLAS 

UN 

-2. 13 

-34.26 

. -0.01 

-36.40 

-36 .47 

-31.14 

7888 

HOPLAS 

UN 

-1.65 

-28.78 

0.00. 

-30.43 

-29.48 

-2 . 20 

7921 

HOPLAS 

UN 

-1.65 

-28.78 

0.00 

-30.43 

-29.48 

-29 .47 

7894 

YUMLAS 

UN 

-3.39 

-30.86 

0.00 

-34.25 

-33.56 

-28 . 63 
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Table 24. Total Geoid Undulations of the Laser Stations Using All the 
Four Methods and the Local Average Correction (Units are in 
meters). 


NUM 

NAME 

STOKES 

MEISSL 

MOLOD. 

SOOB. 

REF 

OSU86F 

ELL-ORT 




WEST UNITED STATES 




705 1 

GUI LAS 

-24.79 

-24.55 

-24.65 

-24.60 

BM 

-23.64 

-22.57 

7109 

QU 1 LAS 

-24.79 

-24.55 

-24.65 

-24.61 

BM 

-23.65 

-22.59 

7886 

GUI LAS 

-24.79 

-24.55 

-24.65 

-24.60 

UN 

-23.65 

-19.30 

7062 

SANOIE 

-34.93 

-34.73 

-34.81 

-34.76 

BM 

-34.22 

-32.47 

71 10 

MONLAS 

-33.02 

-32.70 

-32.84 

-32.78 

BM 

-32.78 

-30.86 

7082 

BEARLK 

-14.36 

-15.11 

-14.78 

-14.57 

BM 

-15.16 

-12.58 

7084 

OVRLAS 

-27.62 

-27.69 

-27.66 

-27.73 

BM 

-26.29 

-24.90 

7114 

OVRLAS 

-27.62 

-27.69 

-27.66 

-27.74 

BM 

-26.30 

-24.84 

7085 

GOLDLS 

-31.79 

-31.74 

-31 .76 

-31.79 

BM 

-31 .70 

-29.68 

71 15 

GOLLAS 

-32.01 

-31.90 

-31.95 

-31.96 

BM 

-31.84 

-29.58 

7086 

FTOAVS 

-20.95 

-21 .06 

-21.01 

-20.96 

BM 

-21.36 

-20.78 

7885 

FORLAS 

-20.95 

-21.06 

-21.01 

-20.96 

BM 

-21 .36 

-20.83 

7112 

PLALAS 

-19.01 

-18.80 

-18.90 

-18.99 

BM 

-18.28 

-17.40 

7887 

VAN LAS 

-36.26 

-36.25 

-36.26 

-36.24 

UN 

-36.47 

-31.14 

7888 

HOPLAS 

-29.99 

-29.84 

-29.90 

-29.94 

UN 

-29.48 

-2.20 

7921 

HOPLAS 

-29.99 

-29.84 

-29.90 

-29.94 

UN 

-29.48 

-29.47 

7894 

YUMLAS 

-33.44 

-33.43 

-33.43 

-33.40 

UN 

-33.56 

-28.63 




EAST UNITED STATES 




706 3 

STALAS 

-32.71 

-32.78 

-32.75 

-32.71 

BM 

-33.08 

-32.92 

706 4 

GSFCLS 

-32.71 

-32.78 

-32.75 

-32.70 

BM 

-33.08 

-33.64 

7100 

GSF 100 

-32.73 

-32.80 

-32.77 

-32.73 

BM 

-33.09 

-32.50 

7101 

GSF101 

-32.73 

-32.80 

-32.77 

-32.72 

BM 

-33.09 

-33.08 

7102 

GSF 102 

-32.71 

-32.78 

-32.75 

-32.70 

BM 

-33.08 

-32.89 

7103 

GSF 1 03 

-32.71 

-32.78 

-32.75 

-32.70 

BM 

-33.08 

-32.93 

7104 

GSF 104 

-32.73 

-32.79 

-32.76 

-32.72 

BM 

-33.08 

-32.94 

7105 

GSF 105 

-32.71 

-32.78 

-32.75 

-32.71 

BM 

-33.08 

-32.93 

706 9 

RAMLAS 

-29.90 

-29.85 

-29.87 

-29.86 

BM 

-29.52 

-29. 13 

709 1 

HAYLAS 

-28.54 

-28.48 

-28.51 

-28.55 

BM 

-28.61 

-27.55 




HAWAII 






7120 

MU I LAS 

18.56 

17.73 

18.10 

18.34 

BM 

13.02 

20.46 

7210 

MAULAS 

18.59 

17.76 

18. 13 

18.37 

IA 

13.02 

20.49 




AUSTRALIA 





7090 

YARLAS 

-24.99 

-25. 18 

-25.00 

-25.05 

BM 

-25.02 

-24.28 

7943 

ORR LAS 

19.39 

19.56 

19.49 

19.48 

UN 

20.36 

20.38 




BERMUDA 





7067 

BD I LAS 

-32.13 

-32. 14 

-32.11 

-32.03 

BM 

-35.23 

-31.81 




EUROPE 






7805 

FINLAS 

20.59 

20.70 

20.64 

20.96 

IA 

20,54 

19.95 

7939 

MATLAS 

46.56 

46.38 

46.46 

45.77 

UN 

45.66 

46.39 

7835 

GRALAS 

51.83 

51.78 

51.80 

51 .47 

UN 

51 .04 

42.85 

78 40 

RGOLAS 

46.31 

46.15 

46.21 

46.39 

UN 

46.06 

45.58 

8833 

KOOLAS 

44.74 

44.78 

44.75 

44.69 

UN 

44.80 

44.66 

7834 

VETLAS 

46.53 

46.87 

46.72 

46.26 

IA 

47.02 

47.65 




JAPAN 






7838 

SHOLAS 

39.56 

39.48 

39.51 

39.11 

UN 

37.60 

39.97 
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Table 25. Differences Between the Undulations Computed Using All the 
Four Methods (and the OSU86F set) and the N T Value (Units 
are in meters). Number of Stations: 28. 


NUM 

NAME 

STOKES 

ME I SSL 

MOLOD. 



OSU86F 

7J0T5 1 

QUILAS 

-1 .78 

-1.77 

-1.77 

-1 .69 

BM 

-1.07 

7109 

QUILAS 

-1.76 

-1.75 

-1.75 

-1.68 

BM 

-1 .06 

7062 

SAND I E 

-2.83 

-2.43 

-2.60 

-2.57 

BM 

-1.75 

7110 

MONLAS 

-2.53 

-2.01 

-2.24 

-2.20 

BM 

-1.92 

7032 

BEARLK 

-4.41 

-3.78 

-4.06 

-4.02 

BM 

-2.58 

7034 

OVRLAS 

-2.52 

-2.69 

-2.62 

-2.67 

BM 

-1.39 

7114 

OVRLAS 

-2.58 

-2.75 

-2 .68 

-2.74 

BM 

-1.46 

7085 

GOIDLS 

-2.26 

-2.13 

-2.19 

-2.23 

BM 

-2.02 

7115 

GOLLAS 

-2.58 

-2.39 

-2.48 

-2.50 

BM 

-2.26 

703 6 

FTDAVS 

-1.39 

-0.86 

-1 .09 

-1.12 

BM 

-0.5 8 

7835 

FORLAS 

-1 .34 

-0.81 

-1 .04 

-1.07 

BM 

-0.5 3 

7112 

PLALAS 

-1 .56 

-1.37 

-1 .46 

-1.55 

BM 

-0.88 

7063 

STALAS 

-0.31 

-0.11 

-0.20 

-0.19 

BM 

-0.16 

7064 

GSFCLS 

0.41 

0.61 

0.52 

0.54 

BM 

0. 56 

7100 

GSF 100 

-0.75 

-0.55 

-0.6 4 

-0.63 

BM 

-0.5 9 

7101 

GSF101 

-0. 17 

0.03 

-0.06 

-0.04 

BM 

-0.01 

7102 

GSF 102 

-0.34 

-0.14 

-0.23 

-0.21 

BM 

-0.19 

7103 

GSF 103 

-0.30 

-0. 10 

-0.19 

-0. 17 

BM 

-0.15 

7104 

GSF 104 

-0.31 

-0. 10 

-0..19 

-0.18 

BM 

-0.14 

7105 

GSF 1 05 

-0.30 

-0. 10 

-0.19 

-0.18 

BM 

-0. 1 5 

7069 

RAMLAS 

-0.09 

-0.40 

-0.26 

-0.20 

BM 

-0. 39 

7091 

HAYLAS 

-1.47 

-1.16 

-1.30 

-1.37 

BM 

-1.06 

7120 

MU I LAS 

-1.53 

-2.55 

-2. 10 

-1 . 83 

BM 

-7.44 

7210 

KAULAS 

-1.53 

-2.55 

-2. 10 

-1.83 

IA 

-7.47 

7090 

YARLAS 

-0.49 

-0.30 

-0.56 

-0.60 

BM 

-0. 74 

706 7 

BD I LAS 

-0.24 

-0.29 

-0.27 

-0.16 

BM 

-3.42 

7305 

FINLAS 

0.53 

0.72 

0.65 

0.96 

IA 

0.59 

7834 

VETLAS 

-0.97 

-0.71 

-0.33 

-1.23 

IA 

-0.63 


Table 26. Differences Between the Undulations Computed Using All the 
Four Methods (and the OSU86F set) with the Local Average 
Correction and the Nj Value (Units are in meters). Number of 
Stations: 28. 


NUM 

NAME 

STOKES 

ME I SSL 

MOLOD. 

SOOB. 

REF 

0SU86F 

7051 

QUILAS 

-2.22 

-1.98 

-2.08 

-2 .03 

BM 

-1 .07 

7109 

QUILAS 

-2.20 

-1.96 

-2.06 

-2.02 

BM 

-1.06 

7062 

SAND I E 

-2.46 

-2.26 

-2.34 

-2.29 

BM 

-1.75 

71 10 

MONLAS 

-2. 16 

-1.84 

-1 .98 

-1.92 

BM 

-1.92 

7082 

BEARLK 

-1.78 

-2.53 

-2.20 

-1.99 

BM 

-2 . 58 

7084 

OVRLAS 

-2.72 

-2.79 

-2.76 

-2.83 

BM 

-1.39 

7114 

OVRLAS 

-2.78 

-2.85 

-2.82 

-2.90 

BM 

— 1.46 

7085 

GOLDLS 

-2.11 

-2.06 

-2.08 

-2.11 

BM 

-2.02 

71 15 

GOLLAS 

-2 . 43 

-2 . 32 

-2.37 

-2.38 

BM 

-2 . 26 

7036 

FTDAVS 

-0. 17 

-0.23 

-0.23 

-0.18 

BM 

-0 . 5 8 

7885 

FORLAS 

-0. 12 

-0.23 

-0. 18 

-0.13 

BM 

-0.5 3 

7112 

PLALAS 

-1.61 

-1.40 

-1.50 

-1 .59 

BM 

-0.88 

7063 

STALAS 

0.21 

0. 14 

0. 17 

0.21 

BM 

—0 . 1 6 

7064 

GSFCLS 

0.93 

0.86 

0.89 

0.94 

BM 

0.56 

7100 

GSF 100 

-0.23 

-0.30 

-0.27 

-0.23 

BM 

-0.59 

7101 

GSF 101 

0.35 

0.28 

0.31 

0.36 

BM 

-0.01 

7102 

GSF 102 

0. 18 

0. 1 1 

0.14 

0.19 

BM 

-0.19 

7103 

GSF 103 

0.22 

0.15 

0.18 

0.23 

BM 

-0.15 

7104 

GSF 104 

0.21 

0. 15 

0. 18 

0.22 

BM 

-0.14 

7105 

GSF 1 05 

0.22 

0. 15 

0. 18 

0.22 

BM 

-0.15 

7069 

RAMLAS 

-0.77 

-0.72 

-0.74 

-0.73 

BM 

-0.39 

7091 

HAYLAS 

-0.99 

-0.93 

-0.96 

-1.00 

BM 

-1 .06 

7120 

MU I LAS 

-1.90 

-2.73 

-2.36 

-2.12 

BM 

— 7.44 

7210 

MAULAS 

-1.90 

-2 . 73 

-2.36 

-2.12 

IA 

-7 . 47 

7090 

YARLAS 

-0.71 

-0.90 

-0.72 

-0.77 

BM 

-0.74 

7067 

BD I LAS 

-0.32 

-0.33 

-0.33 

-0.22 

BM 

-3 . 42 

7805 

FINLAS 

0.64 

0.75 

0.69 

1 .01 

IA 

0.59 

7834 

WETLAS 

-1 . 12 

-0.78 

-0.93 

-1 .39 

IA 

-0.63 
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Table 27. Statistics of the Differences Between the Terrain-corrected N TC 
and Uncorrected N un Undulations and the N T Value Using the 
Modified Sjoberg’s Method for Stations in the Western United 
States (Units are in meters). Number of Stations: 12. 


STATISTICS 

mmssm 

snfsnHi 

Mean Difference 

-2.17 

-2.68 

RMS Difference 

2.31 

2.83 

Standard Deviation Difference 

0.79 

0.92 

Minimum Difference 

-4.02 

-4.32 

Maximum Difference 

-1.07 

-1.29 


Table 28. Statistics of the Differences Between the Undulations Computed 
Using All the Four Methods Nj (and the OSU86F Field) and the 
Nt Value (Units are in meters). Number of Stations: 28. 


STATISTICS 

wmm 

HPB3B 

KJSH 

■3S3B 

Kfisni 


Mean Difference 

-1.27 

-1.18 

-1.21 

-1.21 

-1.39 


RMS Difference 

1.70 

1.65 

1.66 

1.65 

2.36 


Standard Dev. Diff. 

1.13 

1.15 

1.14 

1.12 

1.91 

1.10 

Minimum Difference 

-4.41 

-3.78 

-4.06 

-4.02 

-7.47 

-3.42 

Maximum Difference 

0.58 

0.72 

0.65 

0.96 

0.59 

0.59 


♦excluding stations 7120, 7210 


Table 29. Statistics of the Differences Between the Undulation Computed 
Using All the Four Methods Nj Including the Local Average 
Correction and the N T Value (Units are in meters). Number of 
Stations: 28. 


STATISTICS 

■sn 


WfiS 

HdffiZH 

Mean Difference 

-1.00 

-1.06 

-1.03 

-1.03 

RMS Difference 

1.52 

1.58 

1.54 

1.53 

Standard Deviation Difference 

1.14 

1.17 

1.15 

1.17 

Minimum Difference 

-2.78 

-2.85 

-2.82 

-2.90 

Maximum Difference 

0.93 

0.86 

0.89 

1.01 


8.3 Discussion of the results 

Comparing the terrain-corrected vs. the uncorrected undulations for 
12 stations in the western United States (Table 27) we see that the 
terrain-corrected undulations have better agreement with the N T value 
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than the uncorrected undulations: The mean and the RMS difference 
with the N t value reduced by "50 cm when terrain-corrected data was 
used. The standard deviation of the difference also decreased by "20 
cm when the terrain corrections were taken into account. When we 
compare the four methods with the Nt value for all the 28 stations for 
which we have known height references the standard deviation of the 
difference is "1.10 m (Table 28). From the same Table we can see that 
the standard deviation of the difference between the undulations 
computed using the OSU86F set and the N T values is 1.91 m. If we are 
to exclude the two laser stations 7120, 7210 in Maui that show large 
discrepancies, then the standard deviation of the difference drops to 
1.10 m. If we examine Table 25 we see that there are systematic 
differences for the undulations of the 12 Western United States stations 
and the N T value on the order of 2 m. To further investigate these 
systematic differences the 15'xl5' mean anomalies in the northern 
United States and Canada (Lachapelle et al., 1982) were used to compute 
geoid undulations for the above stations. The statistics of the residuals 
between the values of the undulations using the above data set and the 
N t value are shown in Table 30 for all four methods. For the same 12 
stations the corresponding statistics using the 2'x2' anomalies (Table 
12) are shown in Table 31. We can see that the systematic differences 
are also present when the 15'xl5' anomalies are used. However, the 
standard deviation of the difference decreased from 79 cm (2 x2' data) 
to 65 cm (15'xl5' data) which is surprising since one would expect 
better accuracy when the 2'x2' data are used instead of the 15 x15' 
data. The explanation may be that the 15' averaging procedure was 
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Table 30. Statistics of the Differences Between the Undulations Computed 
Using the 15'xl5' Anomalies and the N T Value for the 12 
Western U.S. Laser Stations. The Modified Sjoberg’s Method 
Has Been Used (Units are in meters). 


STATISTICS 

■SB 

KJBJB 

WEmam 

■UfifBI 

Mean Difference 

-1.99 

-1.72 

-1.83 

-1.82 

RMS Difference 

2.10 

1.83 

1.93 

1.95 

Standard Deviation Difference 

0.65 

0.65 

0.61 

0.72 

Minimum Difference 

-2.88 

-2.66 

-2.64 

-2.69 

Maximum Difference 

-1.04 

-0.73 

-0.87 

-0.62 


Table 31. Statistics of the Differences Between the Undulations 
Computation Using the 2'x2' Anomalies (and the OSU86F Field) 
and the Value for the 12 Western U.S. Laser Stations. The 
Modified Sjoberg’s Method Has Been Used (Units are in 
meters). 


STATISTICS 

BG3I 

bbh 

wmemsm 

WEmm 

■asm 

Mean Difference 

-2.29 

-2.06 

-2.17 

-2.17 

-1.46 

IMS Difference 

2.43 

2.21 

2.31 

2.31 

1.59 

Standard Dev. Diff. 

0.81 

0.81 

0.79 

0.79 

0.63 

Minimum Difference 

-4.41 

-3.78 

-4.06 

-4.02 

-0.53 

Maximum Difference 

-1.34 

-0.81 

-1.04 

-1.07 

-2.58 


better than the 2' averaging procedure. The undulations computed from 
the OSU86F set for the above 12 stations give statistics similar to the 
undulations computed from the 15'xl5' data when compared to the N T 
value. This is an expected result, since the 15'xl5' data was used for 
the computation of the OSU86F set. The cause of the above systematic 
differences was not found and further studies are needed to investigate 
why these systematic differences are present. 

The application of the local average correction improves the 
statistics as it can be seen by comparing Tables 28 and 29. Note the 
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improvement (compare Tables 25 and 26) of about 2 meters of the 
undulation of station 7082 in Bear Lake, Utah when the local average 
correction was used with respect to the N T value. 

Furthermore, there is an excellent internal agreement between the 
undulations computed using all the four methods. The standard 
deviations of the differences between the undulations computed using 
two different methods is on the order of 20 cm. When the undulations 
computed using all the four methods are compared to the undulations 
computed using the OSU86F set they have a standard deviation 
difference of about 1.5 m. The largest differences between the various 
methods and the OSU86F set occur for the 2 laser stations on Maui 
(~5 m) and for the Bermuda station ("3 m). The gravimetric undulations 
for the above cases are closer to the values than the undulations 
computed from the OSU86F set. 

Using a set of 14070 point free-air gravity anomlies by Watts 
(private communication, 1987) in the vicinity of Hawaii the geoid 
undulation of the laser station 7210 was computed to be as follows: 
18.64 m, 17.95 m, 18.26 m and 18.50 m for the Stokes’, Meissl’s, 
Molodenskii’s and modified Sjoberg’s method respectively. We can see 
that the above results agree within ^20 cm with the gravimetric values 
of the undulations of Tables 21 and 24. From the gravity data sent by 
Watts 8.43 mgal were subtracted to account for a systematic difference 
that was detected within the 5*x5* computation area in the vicinity of 
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Hawaii between the Watt’s set and the anomalies implied by the OSU86F 
set. 


8.4 Altimeter Geoid Undulation Computations in Bermuda Area 

An alternative method to compute geoid undulations is the method 
of the least squares collocation. This method was applied to the 10360 
Ohio State adjusted GEOS-3/SEASAT altimeter sea surface heights (Liang, 
1983) covering the Bermuda region (33*<4<31 # .63, 294*.5^X<295'.87, laser 
station number 7067). The above altimeter data was referred to an 
equatorial radius of 6378137 in. To refer the altimeter data to the 

equatorial radius of our adopted reference ellipsoid (see Chapter 4, 
equation (4.1)) we added 1 meter to all the sea surface heights that 
were selected for the collocation method. We then used the collocation 
method to compute the undulation at the laser station 7067 with the 
following specifications: 

a) The whole area (1°.37 x 1".37) was subdivided into 49 0°.20 x 
0*.20 sub-areas so that a maximum of 300 points were used for 
each area by the collocation method. For the predictions near 
the edges of each 0°.20 x 0*.20 area, a border of 0*,20 was 
specified. 

b) The covariance function used between degrees 2 and 180 was the 
one implied by the anomaly error degree variances of the 
potential coefficient solution of Rapp (1981b); between degrees 
180 to ® the Tscherning-Rapp (1974) model (see eq. (2.18)) was 
used for the anomaly degree variances. 
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c) The reference field up to 180 was removed from the data using 
the solution of Rapp (1981b). 

d) The predictions of the undulations were done at the intersections 
of the points of a grid of 2'x2'. The value of the undulation for 
the laser station 7067 was then computed using the four closest 
grid points by bi-linear interpolation with respect to the latitude 
and the longitude. 

The value of the undulation computed from the above collocation 
method was -32.11 m. The accuracy of the undulation predictions at the 
four closest grid points was *0.28 m. 

Using the results of the collocation method and the modified 
Sjoberg’s method two geoid maps (Figures 14 and 15) were produced in 
the Bermuda region (31*.63<4<33*, 294*.5<X*295*.87). The corresponding 
map of the differences between the gravimetric and the altimetric geoid 
(Figure 14 - Figure 15) is shown in Figure 16. The mean difference of 
the 1764 undulations computed at the 2' grid intersections in the above 
region minus the corresponding altimetric undulations was -0.63 m; the 
RMS difference was 0.78 m and the standard deviation of the difference 
was 0.47 m. 

The laser station 7067 has been also used for the calibration of the 
altimeter satellites (Kolenkiewicz et al., 1982). Using the four laser 
supported calibration passes of SEASAT at the time of closest approach 
to the laser station the geoid undulation has been computed for these 
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Figure 14. Gravimetric Geoid Using Modified Sjoberg's Method for the 
Laser Station 7067 in the Bermuda Area (C.I. s 25 cm); 
2'x2' Grid Computed. 


LATITUDE 


dl . b 


31.63 


294.50 295 295.50 

LONGITUDE 


Figure 15. Geoid Computed Using GEOS-3/S EAS AT Altimeter Data for the 
Laser Station 7067 in the Bermuda Area (C.I. = 25 cm); 2'x2' 
Grid Computed. 
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Figure 16. Map of the Differences Between the Gravimetric and the 
Altimeter Geoid in the Bermuda Area (Figure 14 - Figure 15); 
C.I. = 25 cm. 
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passes by Kolenkiewicz et al. (1987), as follows: To the raw altimeter 

data corrections due to tides, dry troposphere, wet troposphere, 
ionosphere, center of gravity, sea state bias and height acceleration 
were applied. Then for the above four closest passes the geoid 

undulation was obtained from the corrected altimetric observations and 
the tracking data by smoothing the offshore altimeter data across land 
(ibid., Table 5). Table 32 summarizes the results obtained from the 
above SEASAT calibration for the four passes: The first columm gives 

the date of pass, the second column gives the latitude, the third column 
gives the geoid undulation which refers to an equatorial radius of 
6378137 m and the last column gives the geoid undulation which refers 
to our adopted equatorial radius of 6378136 m and is obtained by adding 
1 meter to column three. 

Table 32. Geoid Undulations in meters for the Pour SEASAT Passes at 
the Time of Closest Approach to the Laser Station 7067 in 
Bermuda (Kolenkiewicz et al., 1987). 


Date of Pass 
1978 

Latitude 

Undulation 
(a„=6378137 m) 

Undulation 
(a„=6378136 ra) 

Sept. 13 

32 *.3520 

-32.89 

-31.89 

Sept. 16 

32 *.3553 

-32.84 

-31.84 

Sept . 22 

32*. 3615 

-32.69 

-31.69 

Oct. 1 

32 *.3653 

-32.62 

-31.62 


The accuracy that is associated with the above computations (ibid., p. 
27) is *12 cm. To obtain a unique value for the undulation of the laser 
station (♦ = 32*. 3538) we can linearly interpolate this value with respect 
to the latitude from the first two values of the undulation (a e = 6378136 
m) of Table 32. This interpolated value of the undulation at the laser 
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station is then -31.86 m. 

Summarizing the geoid undulation computations for the laser station 
7067 in Bermuda we have: The gravimetric undulations using all the 
four methods (Table 21) are: -32.05 m, -32.10 m, -32.08 m, and -31.97 m 
for the Stokes', Meissl’s, Molodenskii's and modified Sjoberg's methods 
respectively; the "ellipsoidal minus orthometric" value of the undulation 
is -31.81 m; the value computed from the Ohio State GE0S-3/SEASAT 
altimeter data is -32.11 m; and the value resulted from the SEASAT 
calibration is -31.86 m. Thus the results for the Bermuda laser station 
undulation show a very good agreement on the order of 30 cm using all 
the different methods. This is clearly due to the sufficient coverage of 
terrestrial gravity anomalies and GE0S-3/SEASAT altimeter sea surface 
heights in the Bermuda region. 

Finally note that the deviation of the sea surface from the geoid, 
known as Sea Surface Topography (SST) plays an important role in the 
above computations and is separately discussed in Appendix B. 

8.5 Altimeter Geoid Undulation Computations for the Five Laser Stations 
in the Western Atlantic and Pacific Using Collocation 

Finally, the collocation method applied to the adjusted 
GEOS-3/SEASAT surface heights as before was used to compute geoid 
undulations on islands in the Western Atlantic and the Pacific since no 
terrestrial gravity data was available for the above stations. For each 


laser station a window of A$ 


0*.25 in latitude and AX 


0*.25 in 
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longitude centered at the laser station with a border of 0*.25 was used. 
Again the maximum of data points used was 300. The predictions of the 
undulations were done at the intersections of the points of a grid of 
3'x3'. The value of the undulation was then computed using the four 
closest grid points by bi-linear interpolation with respect to the latitude 
and the longitude. Figure 18 shows the undulation map derived from 54 
adjusted GEOS-3/SEASAT altimeter data (whose distribution is shown in 
Figure 17) for the laser station 7061 on Easter Island, Chile. Table 33 
shows the laser station number; the number of adjusted GE0S-3/SEASAT 
altimeter data that were used in the predictions for each oceanic station; 
the undulation N ALT computed from the OSU86F set taken up to degree 
360; and the differences N alt -Nt and N P0T -N T . 

In Table 34 we can see the statistics of the differences N AL t~Nt and 
NpofNi f°r the five oceanic stations. The results using the collocation 
method showed better agreement with the N T value than the OSU86F set. 
The standard deviation of the difference in the first case is 0.79 m and 
in the second case is 1.25 m. Also the mean difference in the first case 
is by 1.85 m smaller than in the second case. The better agreements 
with the N t values using the collocation method than using the OSU86F 
set can be explained as follows: Within the collocation method used here 
all the available altimeter sea surface heights in a 0*.75 x 0*.75 area 
were taken into account, whereas for the OSU86F solution 0*.5 x 0*.5 
mean anomalies derived form the above sea surface heights were used, 
where clearly some loss of the high frequency information has taken 
place. 
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Figure 17. Distribution of the 54 Adjusted GEOS-3/SEASAT Altimeter 
Data Surrounding the Laser Station 7061 in Easter Island, 
Chile U+=AA=0*.75). 
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Figure 18. Geoid Undulation Map Based on the Altimeter Data of Fig. 
17. C.I. = 5 cm (A4>=AA=0*.20). 
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Finally, the effect of the Sea Surface Topography is separately 


discussed in Appendix B. 


Table 33. Geoid Undulations Computed Using the Collocation Method 
(N AL t)> Based on Adjusted GEOS-3/SEASAT Altimeter Data and 
from the OSU86F Field, for the Five Oceanic Laser Stations 
(Units are in meters). 


STATION 

DATA 

n all 

NpoT 

- Nt - 

BMVUifli 

m&nii 


54 

- 3.26 

- 3.35 

- 2.25 

-1.01 

-1.10 

■ 

299 

-42.63 

-44.81 

-42.53 

-0.10 

-2.28 

■ 

127 

31.77 

29.72 

32.18 

-0.41 

-2.46 

7096 

40 

33.36 

30.45 

34.69 

-1.33 

-4.24 

7121 

96 

8.00 

5.96 

10.35 

-2.35 

-4.39 


Table 34. Statistics of the Differences N ALT -N T and N PO t - Nt f° r the Five 
Oceanic Stations (Units are in meters). 


STATISTICS 

N AL t-Nt . 

Npot-Nt 

Mean Difference 

-1.04 

-2.89 

RMS Difference 

1.30 

2.04 

Standard Deviation Difference 

0.79 

1.25 

Minimum Difference 

-2.35 

-4.39 

Maximum Difference 

-0.10 

-1.10 









CHAPTER IX 


SUMMARY AND CONCLUSIONS 

The geoid undulations for 39 laser stations distributed around the 
world have been computed using the classical Stokes’ and Meissl's 
methods and the new methods of Molodenskii’s and modified Sjoberg’s. 
The high frequency information for the undulation was taken from 

gravity data within a cap of radius 2* and the low frequency 
information was taken from a set of potential coefficients (OSU86P). 
Using the ellipsoidal heights of the laser stations available from the SL6 
satellite dynamic solution and the orthometric heights H, available from 
spirit leveling, the values of the geoid undulation N T =h-H was obtained 
and served as a mean for comparisons. 

The theoretical improvement in the accuracy of the geoid undulation 
computations using the modified Sjoberg’s method was on the order of 
10 cm (Figure 1) which is well below the noise of the actual data used 
and thus this theoretical improvement was not realized in practice. No 

method was clearly better than the others. However, all the four 

methods gave internally consistent results on the order of 30 cm and 
for most of the regions (Australia, Europe, Bermuda, eastern United 
States) there was an agreement of ~60 cm or better between the 

undulations computed from the four methods and the N T value (see Table 
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25). For some regions where local high variation of the gravity field 
occurs (e.g. Maui, (stations 7120,7210), Bermuda (station 7067)) the 
OSU86F field gives high discrepancies as compared to the N T values. 
The "worst" situation was in Maui, where the OSU86F value of the 
undulation differs from the N T value by about 7 meters (Table 25). This 
large difference can be explained by the fact that the OSU86F, being a 
limited expansion to degree 360 cannot represent the very high 
frequency signal of some regions of high variation of the gravity field. 
Such regions are the islands of Maui and Bermuda. The effect of the 
terrain was taken into account for the 17 laser stations in the Western 
United States where the terrain was considered to be rough, The effect 
of this correction was on the order of 60 cm (Table 27). For most of 
the laser stations no terrain corrections and/or height information was 
available and thus the computed value of the undulation will be 
systematically low. This is clearly the situation for the two stations in 
Hawaii, where the difference between the computed undulations using 
the four methods and the N T value is -2 m (see Table 25). The 
ellipsoidal corrections for the four methods were derived in a uniform 
formulation. Using the OSU86F set up to degree 36 those corrections 
were computed and they were very small in magnitude (maximum: 3 cm 
which is in good agreement with the results obtained in (Rapp, 1981b, 
Table 1)). Consequently, they could be neglected in the total undulation 
computation. A very important correction due to inconsistencies between 
the local terrestrial anomalies and the anomalies implied from the 
geopotential set within the same region of a 2* cap, called the "local 
average correction" was applied to all the computed undulations. The 
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effect of this correction reached 2 m for station 7082 in Bear Lake, Utah 
(compare Tables 25,26) and it improved the results with respect to the 
N t value by about 30 cm for 28 laser stations that were used for the 
comparisons (Tables 28,29). It seems that this correction has to always 
be applied when geoid undulations are computed by a combination of 
terrestrial data and potential coefficient information. Very good 

agreement (**30 cm) was achieved between the gravimetric undulations 
using the four methods, the altimetric undulation, the undulation 
obtained from the calibration of SEASAT (Kolenkiewicz, 1987) and the N T 

value for the station 7067 in Bermuda (see discussion at the end of 

Section 8.4). 

The collocation method using Ohio State adjusted GE0S-3/SEASAT 
altimeter data for five oceanic stations resulted in values of the 
undulation consistent with the N T value on the order of 70 cm (Table 
34). 

Using the low of propagation of errors the accuracy of the 

computed geoid undulations for each of the 44 laser stations could be 

computed from the errors in the gravity data surrounding the laser 
station, the accuracy of the given potential coefficients and from various 
models that could represent the omission and discretion errors. These 
accuracies were not computed individually for each laser station 
although global error models were considered in the error analysis of 
the various methods in Chapter 3. The global RMS undulation error 
from Chapter 3 for the capsize of 2* was computed to be approximately 
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60 cm for all the four methods (Figure 1). Finally, the accuracies of the 
altimetric undulations for all the five oceanic laser stations using the 
collocation method were computed to be ~20 cm. 


CHAPTER X 


RECOMMENDATIONS FOR FUTURE STUDIES 

It is very critical for the future applications of the modified 
Sjoberg’s method that error degree variance models of the terrestrial 
data which will represent the actual local behavior of the error 
estimates of those data should be developed. The modification of the 
Stokes’ function for the modified Sjoberg’s method is essentially based 
on the above error models and if the models are properly chosen this 
method will give the smallest possible error in the geoid undulation 
computations. A Monte-Carlo method modeling the empirical error 
covariance function could be suggested. Then, the error degree 
variances can be obtained by integration of the error covariance 
function. 

The anomaly error degree variances implied by a geopotential model 
would certainly be different if a full variance-covariance matrix of the 
potential coefficients could be available. This would again influence the 
modified Sjoberg's method but the exact effect was not tested. 

Another important investigation for precise undulation computations 
can be the effect of the zero and first degree term of the error anomaly 
degree variances of the terrestrial data and the zero degree term of the 
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error degree variances of the potential coefficient data (<r 0 , <r x and 5C 0 
respectively, following the notation of the introduction) on the 
computation of the W 4n coefficients in the modified Sjoberg’s method. 
Although the theory was given by equation (2.106) and (2.107), for the 
numerical applications we assumed a , o = £C o =<ri=0. It seems that the 
stability of the solution (2.108) will be increased by admitting zero and 
first degree term for the above error degree variances, although this 
was not checked. The local average correction should be more 

extensively checked using actual geoid computations. Gravimetric geoid 
undulation computations in Doppler or GPS stations including and 
excluding the zero order correction could give a definite answer 
whether this correction should be applied or not. However, the answer 
reached here is that this correction has to always be applied. 

Finally, the collocation method of computing undulations using the 
randomly distributed point gravity data surrounding the laser stations 
could be applied and the results could be compared with the ones given 
here using the integral formulas. 
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APPENDIX A 


Proof of equations (2.105), (2.106) and (2.107): 

The expression of the global mean square error combining equations 
(2.98), (2.99), (2.100), (2.101) and (2.102) is: 

+ ^ , (n^l ~ Q *"Y C " + 8 (Q-n + W 4n ) 2 5C n 

n = n T +i '• n 1 1 n = o 


+ X. Q5n«G n + X <tf„c n ] 

n-n+i n-Mfi - 


/ K 1 \ 

V/i. J.y 


For convenience, let us use the substitutions 


X«n = X n 

W*n = W n 


(A. 2) 


_R 

2y 


c 


(A. 3) 


The minimization of 6N| expressed in (2.103) is equivalent to the 


conditions 


. o, m.i. — s 


(4.4) 


where, from (A. 1) and (A, 2), (A. 3): 

«S = c*[ l (Xn-a.) 1 ’. + V. (^rx - + i (jrr 

L n=o n=n+i vn 1 } n=n T +i vn 1 

+ X (Q n + W„) a 6C n + l Q 2 fiC n + i Q 2 C n ] 

n = o « = f.rU-4 - 1 J 


~ Qn ) 2 C n 


n = n+i 


n-H+i 


(A. 5) 


Denoting with Xi» X2, Xs> X4, Is and X« the summations in (A.5) and 
taking into account eq. (2.96), we have: 
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h = I (X n -Qn) 2 ^n = X (x„-Q ln + X ^ W r e rn ) 2 ff n 
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Thus, condition (A. 4) can be written 


3(g W2) _ 3fc a (Z,+i: a +I,+I^E,+j 
dw k aw k 

_ c 2 3(y,+y, +£,+£. +v,+y ( 
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or, since c 2 *Q 
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3(Q„) _ _ 2k+l 

3W k ~ 2 ekn » k= °.l» . . . ,n 

where fi kn is the Kroenecker’s 6: fi kn = 
partials (A. 7) from (A. 6), one by one: 


we evaluate the 
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where, for (A. 9 ) - (A. 14 ) K= 0 , 1 , . . . ,ii. 
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the partials (A. 9) - (A. 14) can be written: 
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(A. 17 ) 



124 


, 2k+l St S 2r+l B 

9 ? 2 ®r e rn e kn^n 

1 2k+l “ _2_ 2k+l ST . _ 

2 9W k - 2 n= i T+1 n-1 S| <" C " 2 „h+V” ek ” C ” 

2k+l • n 

+ 1 , 1 W rernCn 

^ n=nj+l r-0 

2 “ ~~ 2 ~ ? ®i n e kn^n + E Qin^kn^n 

^ an k ^ n=o n=o 

. 2k+l 9 9 2r+l M 

2 ^ 2 ^r e rn e kn^n 

^ n=0 r=0 ^ 

n n 2r+1 n n 

~~ ? ? ^r e rn^kn^n “ T" E e knW n ^kn + E W n <5 kn 6C n 

r-o r-o ^ ^ n =0 n=o 

1 l£i - _ 2k+l H 2k+l M n 

2 3Wu 2 ® 1 n e kn^C n + g X- 2 ^r e rn e kn^n 

6 n k ^ n=n+i ^ n— n+l r=o 

1 il* _ _ 2k+i 9 2k+i 9 M 

2 3W k ■ 2 n =5 +1 QlneknCn 2 n= L 1 Wrer " ek " C " 

Using condition (A. 7) and adding (A. 17) - (A. 22), we obtain 
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? n — 1 ^kn®n + ? ®in®kn^n + I Qln^kn ff n ~ 0, k-0,1, . . . ,n(A.23) 

n — 2 11 * n — o n=o 

The third term of the left hand side of eq. (A. 23) can be written: 

n n 9r+l n o -.- 1 i 

^X^ 2~ ^r e rn^kn( a n + ^n) = X^ 2 ^r e rk(^k + ^k) (A. 24) 

Using (A.24) into (A.23) we obtain after collecting the common factors of 
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where X lk are given by (2.39). 


It is now obvious that equations (A.25) are equivalent to the linear 
system 

X a kt^t ” ^k k=0, 1, . . . ,n (A. 26) 
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+ n=M+» ekn [® ln ^ ff " + ^ n ^ ~ 
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n-ny+i 1 JJ 

+ X lk <r k - Q lk (<r k +5C k ) (A. 28) 

Equations (A.26), (A.27) and (A.28) after the substitutions (A.2) and (A.3) 
give immediately equations (2.105), (2.16) and (2.107) which were to be proven. 


APPENDIX B 


The deviation of the sea surface from the geoid, known as the Sea 
Surface Topography (SST) plays an important role when the gravimetric 
or the altimetric undulation are compared to the ellipsoidal minus 
orthometric value of the undulation of laser stations on islands. From 
Figure B1 we can see that the orthometric height (H) of a laser station. 



Figure Bl. Role of the Sea Surface Topography (£) for Undulation 
Computations in Oceanic Laser Stations. 


the SST (£), the geoid undulation (N') and the ellipsoidal height (h) of 
the laser station satisfy the equation: 
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h = H + f + N' (B.l) 

In writing equation (B.l) the basic assumption that has been made is 
that the vertical datum to which the orthometric height of the laser 
station refers is assumed to coincide with the mean sea surface as is 
realized by the altimetric observations. From (B.l) the value of the 
undulation will be 

N' = (h-H) - f (B. 2) 

From Chapter 4, equation (4.2), and from equation (B.2) we take: 

N' = N T - f (B.3) 

Thus, from the value of the undulation computed in Chapter 4, the 
value of the SST has to be subtracted, so that the undulation N' 
obtained can be compared to the gravimetric undulation of the laser 
station. We thus assume that the gravimetric undulations of the laser 
stations are free from the effect of SST. This is not rigorously true, 
since the SST influences the gravity anomaly data which are used to 
compute the gravimetric undulations, but this effect can be considered 
negligible (Rapp, 1985b). 

Four different SST harmonic coefficient sets described in (Rapp, 
1985b, Appendix B) were used to compute f. These sets were computed 
by Engelis (private communication, 1987) similarly as in (Engelis, 1985b) 
and they are as follows: 

Set No. 1 : SST harmonic coefficients complete to degree and order 36. 

The PGS3041 potential coefficient set, an intermediate set for 
the GEM-T1 (Marsh et al., 1986) has been used to realize the 
geoid up to degree 36. The high frequencies from degrees 
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37 to 180 were removed from the SEASAT sea surface heights 
by using the potential coefficients of the solution described 
in (Rapp, 1981b). 

Set No. 2 ; SST harmonic coefficients complete to degree and order 36; 

the PGS3041 set has been used up to degree 36. 

Set No. 3 ; SST harmonic coefficients complete to degree and order 10, 
based on the Levitus (1982) data. 

Set No. 4 : A set of 64800 l*xl* global values of SST based on Levitus 
data. 

The computation of the SST at the laser stations using Sets No. 1-3 was 
done as follows: 

M n 

f = I I (c nm cosmX + d nm sinmA)P nm (sin$) (B.4) 

n-o m-o 

where c nm , d nm are the SST harmonic coefficients and $ is the geocentric 
latitude of the laser station. The degree of expansion M was taken 6 
for all the three sets. 

The computation of the SST at the laser stations using the Set No. 4 
was done using the four l*xl* closest points to the laser station in a 
bi-linear interpolation with respect to the latitude and longitude of the 
laser station. 

The results of the SST computations using all the four SST sets 
described above are shown in Table B1 for 5 laser stations. The 
number and the name of the laser stations appear in the first two 
columns; the location of the stations appears in the third column; finally 
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the SST as computed from Sets No. 1-4 ((S1),(S2),(S3),(S4) 

correspondingly) is shown in columns four-seven. The units for the 
SST is meters. 


Table Bl. Sea Surface Topography in meters for 5 Laser Stations Using 
the Sets No. 1-4. 


1 NUMBER 

NAME 

OCCUPATION 

SI 

S2 

S3 

S4 


MAULAS 

Maui, Hawaii 

da 


0.44 

0.32 

I 

BDILAS 

Bermuda Island 

mm 

mm 

0.07 

0.09 

7838 

SHOLAS 

Simosato, Japan 


E8H 

0.43 

0.22 

7090 

YARLAS 

Yarragadee, Australia 

-0.10 

-0.12 

■ail 

0.00 

7943 

ORRLAS 

Orroral Valey, Australia 

0.16 

0.10 


0.00 


For the above five laser stations the gravimetric undulations are 
available using four different methods (See Table 21). We select here 
the undulations computed using the modified Sjoberg’s method to 
represent the gravimetric undulation of stations 7210-7943 of Table Bl. 
These undulations are shown on the second column of Table B2‘, then, 
the N t value of the undulation computed from (4.2) is shown in the 
third column; the corrected values of N T through equation (B.3) using 
the above four SST sets (N'=N-|--0 is shown in columns four-seven. 
According to our discussion above the undulations of column two should 
be compatible only with the values of the undulations of columns 
four-seven of Table B2. 


The differences N-N T and N-N' for the four SST sets are shown on 
Table B3; Table 37; the corresponded statistics of these differences are 
shown on Table B4. It can be seen from Table B4 that the mean and 


the RMS differences were decreased in all cases that the SST correction 
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Table B2. Gravimetric Value of the Undulation (N), N T (Uncorrected for 
SST) and N'=Nt (Corrected for SST Using the Sets No. 1-4) 
for 5 Laser Stations in meters. 


NUMBER 

N 

Nt 

KHSn 

1 mti 1— 

mamam 

BB&nt 

7210 

18.66 

20.49 

19.68 

19.61 

20.05 



-31.97 

-31.81 

-32.36 

-32.42 

-31.88 

-31.90 

7838 

38.63 

39.97 

39.04 

39.11 

39.54 

39.75 


-24.88 

-24.28 

-24.18 

-24.16 

-24.39 

-24.28 

I 7943 

18.85 

20.38 

20.22 

20.28 

20.41 

20.38 


was applied to the Nj values as compared to the mean and the RMS 
differences with the uncorrected values. The standard deviation of 
the differences remained essentially the same when the first two SST 
sets were used, and slightly decreased (~0.05 m) when the last, two SST 
sets were used. 


The accuracy of the SST computed using the four sets as above is 
thought to be approximately *20 cm (Rapp, 1985b, p. 106). Considering 
this accuracy estimate is true, we can say that the inclusion of the SST 
effect clearly improved the individual and overall comparisons of the 
gravimetric undulations with the (corrected) N T value. 


Table B3. Differences N-N T , N-N'(Sl), N-N'(S2), N-N'(S3) and N-N'(S4) in 
meters (see Table 36) for the 5 Laser Stations. 


NUMBER 

■SEIM 



ESiffE&a 



-1.83 

-1.02 


-1.39 

-1.51 


-0.16 

0.39 

0.45 

-0.09 

-0.07 

7838 

-1.34 

-0.41 


-0.91 

-1.12 

7090 

-0.60 

-0.70 


-0.49 

-0.60 

7943 

-1.53 

-1.37 

-1.43 

-1.56 

-1.53 
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Table B4. Statistics of the Differences N-N T> N-N'(Sl), N-N'(S2), N-N'(S3) 
and N-N'(S4) in meters (see Table 36) for the 5 Laser 
Stations. 


STATISTICS 

N-N t 

UffijftttM 



IfflUflSEiM 

Mean Difference 

-1.09 

-0.62 

-0.63 

-0.89 

-0.97 

RMS Difference 

1.25 

0.86 

0.88 

1.04 

1.12 

Standard Dev Diff 

0.61 

0.60 

0.61 

0.54 

0.56 

Minimum Difference 

-1.83 

-1.37 

-1.43 

-1.56 

-1.53 

Maximum Difference 

-0.16 

0.39 

0.45 

-0.09 

-0.07 


Although the sample of the 5 laser stations used here is small, the 
systematic improvement of the results when the SST was considered is 
clear. We thus can say that the SST effect has to always be taken into 
account when geoid undulations at an accuracy level of 50 cm or better 


are to be computed. 







